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ABSTRACT 


Objectives :  Predicting  the  environmental  fate  of  novel  munitions  compounds  requires 
data  for  the  key  properties  that  control  contaminant  fate,  and  the  most  common  way  to 
predict  these  properties  is  with  empirical  quantitative  structure-activity  relationships  (e.g., 
QSARs).  However,  the  traditional  approach  to  QSAR  development  (calibration  with 
experimental  data)  is  challenging  for  new  explosives  compounds  because  of  limitations 
to  the  availability  of  these  materials.  Some  of  the  necessary  environmental  fate  properties 
can  be  calculated  directly  from  molecular  structure  theory,  but  reliable  calculations  of 
this  type  require  considerable  theoretical  expertise  and  computational  effort.  To 
overcome  this  combination  of  challenges,  we  used  a  hybrid,  partially  in  silico  approach  to 
QSAR  development,  where  some  of  the  calibration  data  (both  the  target  variable  and  the 
descriptor  variables)  were  calculated  from  molecular  structure  theory. 

Technical  Approach :  The  technical  approach  of  the  project  we  divided  into  three 
objectives,  corresponding  to  the  three  reaction  pathways  that  contribute  most  to  the 
overall  fate  of  explosives  in  the  environment.  These  include:  (i)  hydrolysis  and  associated 
elimination  reactions,  which  are  ubiquitous  in  water;  (ii)  nitro  reduction  by  outer-sphere 
electron  transfer,  which  are  dominant  in  soils  in  sediments;  and  (iii)  oxidation  of  the 
amino  products  soil  minerals  like  manganese  dioxide. 

Results  '.  After  extensive  analysis  of  prior  work,  it  was  concluded  that  mechanism(s)  of 
hydrolysis  of  the  nitro  aromatic  energetic  compounds  TNT  and  DNAN  was  not 
sufficiently  advanced  to  support  conventional  QSAR  development.  The  possible 
difference  in  mechanisms  means  that  predicting  the  reaction  mechanism  for  one  based  on 
the  other  may  lead  to  unreliable  predictions  of  environmental  fate.  This,  along  with 
uncertainties  in  the  consistency  of  the  calculated  results  with  experimental  values, 
presents  a  challenge  for  developing  QSARs  calibrated  “m  silico ”  that  predict  the 
hydrolysis  behavior  of  the  diverse  range  of  energetic  NACs.  However,  new  experimental 
and  computational  results  reported  here  provide  insight  into  these  mechanisms. 

In  contrast  to  the  hydrolysis  pathway,  the  effect  of  complexity  in  the  nitro  reduction 
pathway  was  effectively  captured  by  contrasting  two  approaches:  theory-based 
correlations  to  Ex  and  empirical  correlations  to  Z/lumo-  The  former  provides  a  variety  of 
insights  into  the  fundamental  processes  controlling  nitro  reduction  and  the  latter  provides 
simple  QSARs  that  should  be  useful  for  estimating  relative  reduction  rates  of  nitro 
aromatic  compounds.  The  models  were  calibrated  and/or  tested  with  a  variety  of  model 
and  actual  energetic  compounds,  including  TNT  and  DNAN.  An  overall  conclusion  of 
this  aspect  of  the  study  is  that  most  alternative  energetic  compounds  will  react  by  nitro 
reduction  more  slowly  than  TNT. 

For  oxidation  of  anilines  a  variety  of  QSARs  were  developed  that  represent  a  range  of 
convenience  and  accuracy.  The  QSARs  based  on  the  most  convenient  descriptors  (e.g., 
pKa)  were  less  accurate  and  models  based  on  the  more  specialized  descriptors  (e.g.  ,  ZZ ' ox) 
were  the  most  accurate.  Comparison  of  the  model  for  MnCh  to  other  environmentally- 
relevant  oxidants  showed  that  the  fonner  was  more  sensitive  to  substituent  effects  than 
the  latter. 

Benefits.  The  QSARs  provided  by  this  project  can  be  used  to  predictive  the  relative  rates 
of  major  transformation  processes  that  contribute  to  the  environmental  fate  of  energetic 


compounds.  The  models  are  flexible  enough  to  handle  some  of  the  diverse  molecular 
structures  that  are  found  in  candidate  compounds  for  new  explosives.  They  are  being 
integrated  into  models  developed  by  others  to  predicted  the  overall  fate  of  energetic 
compounds. 
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OBJECTIVES 


Predicting  the  environmental  fate  of  novel  munitions  compounds  requires  data  for  the  key 
properties  that  control  contaminant  fate,  and  the  most  common  way  to  predict  these 
properties  is  with  empirical  quantitative  structure-activity  relationships  (e.g.,  QSARs). 
However,  the  traditional  approach  to  QSAR  development  (calibration  with  experimental 
data)  is  impractical  for  new  explosives  compounds  because  of  limitations  to  the 
availability  of  these  materials.  Some  of  the  necessary  environmental  fate  properties  can 
be  calculated  directly  from  molecular  structure  theory,  but  reliable  calculations  of  this 
type  require  considerable  theoretical  expertise  and  computational  effort. 

To  overcome  this  combination  of  challenges,  we  proposed  a  novel,  in  silico  approach 
to  QSAR  development,  where  most  of  the  calibration  data  (both  the  target  variable  and 
the  descriptor  variables)  are  calculated  from  molecular  structure  theory.  Rates  of  the  most 
likely  breakdown  pathways  (target  variables)  were  to  be  calculated  with  the  highest  level 
of  theoretical  accuracy  and  these  data  were  correlated  to  molecular  properties  (descriptor 
variables)  that  were  obtained  with  computational  methods  that  are  more  available  and 
feasible  for  most  chemists.  As  QSARs  were  obtained  by  this  approach,  we  attempted  to 
validate  them  by  predicting  data  for  (safe  and  available)  model  compounds  and 
comparing  to  measured  experimental  values. 

The  technical  approach  of  the  project  was  originally  divided  into  four  objectives  (later 
streamlined  to  three),  corresponding  directly  to  the  three  reaction  pathways  that 
contribute  most  to  the  overall  fate  of  explosives  in  the  environment.  These  four  pathways 
are: 

•  Hydrolysis  and  associated  elimination  reactions,  which  are  ubiquitous  in  water. 

•  Homogeneous  nitro  reduction  by  outer-sphere  electron  transfer  with  dissolved 
electron  donors  or  at  the  surfaces  of  reducing  minerals. 

•  Polymerization  of  the  amino  products  (with  themselves  and  with  the  phenolic 
moieties  associated  with  natural  organic  matter)  of  reduction,  by  oxidative 
coupling  or  nucleophilic  condensation. 

Each  one  of  these  objectives  was  pursued  via  five  tasks:  (/)  reaction  formulation  and 
mining  of  existing  data,  (ii)  calculation  of  target  variable  data  with  high  level  theory,  (iii) 
calculation  of  descriptor  variable  data  with  lower  level  theory,  (iv)  correlation  analysis 
and  fitting  QSARs,  and  (v)  validation  of  data  predicted  with  the  QSAR. 
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BACKGROUND 


Predicting  chemical  properties  is  a  long-standing  challenge  that  has  received  extensive 
study  for  many  applications  (chemical  engineering,  green  chemistry,  environmental 
chemistry,  toxicology,  pharmacology,  etc.).  For  chemical  properties  that  determine  the 
fate  and  effects  of  environmental  contaminants,  the  field  is  mature  enough  to  have 
already  engendered  several  generations  of  compilations  of  predictive  models 1 .  However, 
a  realistic  assessment  of  the  state-of-the-art  in  this  field  reveals  some  pervasive 
limitations.  In  particular,  most  existing  property  prediction  models  are  based  on  empirical 
correlations  between  the  target  property  and  descriptor  variables  (e.g.,  quantitative 
structure-activity  relationships,  QSARs),  which  must  be  calibrated  with  measured  data, 
and  are  of  dubious  validity  for  extrapolation  outside  the  training  set  data-space. 
Expanding  or  validating  the  scope  of  application  of  such  QSARs  requires  more 
experimental  data,  which  often  ends  up  being  the  same  data  that  the  QSAR  was  supposed 
to  serve  to  estimate.  To  overcome  this  dilemma,  novel  approaches  to  chemical  property 
prediction  are  needed. 

With  respect  to  chemical  property  prediction  to  support  development  of  new, 
environmentally  benign  chemicals  for  explosives,  the  challenge  noted  above  is 
particularly  severe  for  two  reasons.  First,  safety,  regulatory,  and  security  considerations 
make  it  impractical  to  do  the  experiments  necessary  to  calibrate  and  validate  QSARs. 
Second,  while  most  current  explosives  and  candidates  for  future  explosives  are 
characterized  by  >N-NC>2  or  Ar-NCfi  moieties  [/],  they  do  not  fit  into  simple  homologous 
series  of  congeners  that  are  described  by  traditional  QSARs.  This  can  be  seen  in  Figure 
1,  where  we  have  shown  the  structures  of  some  current  and  future  explosives  arranged  to 
illustrate  groups  that  have  a  high  enough  degree  of  homology  and  number  of  potential 
members  to  make  meaningful  QSARs  feasible.  The  groups  are  defined  by  a  variety  of 
analogous  backbone  types  containing  variable  numbers  of  a  single  moiety  that  is  subject 
to  a  common  pathway  of  reaction.  The  properties  of  such  groups  of  chemicals  can  be 
described  by  QSARs,  but  they  must  be  based  on  suitable  descriptor  variables.  Most 
descriptor  variables  that  are  suitable  for  this  purpose  are  molecular  properties,  which  are 
generally  measured,  which — as  noted  above  for  target  properties — is  impractical  for  most 
explosives. 

Alternatively,  most  properties  that  directly  impact  environmental  fate  can — in 
principle — be  calculated  directly  from  molecular  structure  theory.  However,  perfonning 
accurate  calculations  of  this  type  requires  specialized  theoretical  expertise  and 
considerable  computational  resources,  so  the  fully  uab  initio ”  approach  to  chemical 
property  estimation  is  also  not  a  practical  solution  for  most  of  those  who  are  responsible 
for  assessing  the  potential  environmental  fate  and  effects  of  novel  explosives. 

To  overcome  this  unique  combination  of  challenges,  we  proposed  take  a  novel  hybrid 
approach  that  provides  QSAR-like  models  calibrated  with  data  obtained  partially  from 


1  Most  prominent  among  these  is  the  Handbook  of  Chemical  Property  Estimation  Methods  compiled  by 
Lyman  et  al.  [2],  which  became  so  highly  regarded  that  a  subsequent  volume  by  Mackay  and  Boethling  [3] 
was  promoted  as  the  2nd  Edition  of  Lyman  et  al.  and  its  publication  was  consecrated  with  very  well 
attended  symposium,  with  Lyman  as  keynote  speaker,  at  the  20th  Annual  meeting  of  the  Society  for 
Environmental  Toxicology  and  Chemistry,  14-18  November  1999,  Philadelphia,  PA. 
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theoretical  calculations  based  on  molecular  structure  theory.  This  “fully  in  silico ” 
approach  to  QSAR  development  uses  descriptor  variable  data  that  can  be  calculated  using 
methods  that  are  feasible  for  non-specialists  and  target  variable  data  calculated  (instead  of 
measured)  for  complete  characterization  of  reaction  coordinates  using  state-of-the-art 
computational  methods  of  highest  accuracy.  For  example,  consider  the  formation  of 
hydride-Meisenheimer  complexes,  which  frequently  controls  rates  of  biodegradation  for 
2,4,6-trinitrotoluene  (TNT)  and  other  nitro  aromatics  [4],  There  are  no  reliable  kinetic 
data  for  this  reaction  and  they  would  be  challenging  to  obtain  experimentally,  but  rates 
constants  for  this  reaction  could  be  calculated  from  theory  since  the  reduction  mechanism 
is  well-defined.  Correlating  the  calculated  rate  constants  (target  variable)  with  an  easily 
calculated  molecular  property  (descriptor  variable)  like  charge  density  on  key  atoms 
should  produce  a  QSAR  that  would  be  easy  to  apply  to  other  nitro  aromatics. 
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Figure  1.  Structures  of  key  members  of  three  families  of  modern/emerging  energetic 
substances,  where  each  family  represents  a  series  of  compounds  that  are  sufficiently 
closely  related  to  form  the  basis  for  a  QSAR  model  of  an  environmental  property. 


In  principle,  this  approach  could  be  applied  to  any  process  affecting  the  fate  of  any 
type  of  contaminant,  but  the  focus  of  this  study  on  explosive  compounds  allowed  us  to 
narrow  the  range  of  processes  to  be  considered.  By  necessity,  explosive  compounds  are 
generally  oxidized,  polar,  and  unstable  compounds.  Therefore,  in  the  environment,  they 
show  relatively  little  tendency  to  adsorb  or  bioaccumulate,  facile  reduction  but  relatively 


5 


little  susceptibility  to  oxidation,  and  significant  hydrolysis  initiated  by  “spontaneous” 
elimination  reactions.  This  means  that  we  can  cover  the  most  important  processes 
controlling  the  environmental  persistence  of  explosives  by  considering  just  hydrolysis 
and  reduction  reactions.  Hydrolysis  can  be  adequately  described  as  a  reaction  in 
homogeneous  solvent,  but  modeling  reduction  is  more  complex.  Since  the  overall  rate  of 
most  contaminant  reductions  is  controlled  by  the  rate  of  the  initial  electron  attachment, 
we  developed  a  generic  model  for  reduction  by  single  electron  transfer  from  a  solution- 
phase  donor. 

Reduction  of  nitro  groups  produces  amino-substituted  products  and  the  environmental 
fate  of  these  amines  is  controlled  mainly  by  oxidative  coupling  to  soil/sediment  organic 
matter  [5,  6].  Like  reduction,  this  is  a  complex  and  potentially  system-specific  process, 
but  we  should  be  able  to  model  relative  rates  of  the  process  by  assuming  a  rate-limiting 
oxidation  to  produce  radical  intermediates  that  couple  with  phenolic  moieties  of  natural 
organic  matter  (NOM)  model  compounds.  In  this  case,  disappearance  of  the  original 
contaminant  is  not  the  most  challenging  issue;  instead,  the  greatest  concern  is  with 
assessing  the  potential  for  reversing  this  reaction  (decoupling)  and  release  of  the  amino 
residues,  which  are  potentially  problematic  contaminants.  With  the  theoretical  approach 
used  here,  it  should  eventually  be  possible  to  consider  both  coupling  and  decoupling  (i.e., 
the  reversibility  of  this  process). 

While  all  of  the  model  development  and  some  of  the  model  validation  used  in  this 
study  were  done  in  silico,  we  tested  each  model  against  experimental  data  in  several 
ways.  First,  where  experimental  data  were  available  for  the  target  or  descriptor  variables 
describing  hydrolysis,  reduction,  or  oxidative  coupling  of  explosives,  we  compared  these 
directly  with  our  calculated  values.  However,  since  there  were  many  data  gaps,  our 
approach  had  the  fundamental  advantage  that  we  can  calculate  properties  of  other  model 
compounds  for  which  data  are  available  and  use  these  in  validation  of  our  QSARs. 
Finally,  we  used  the  combination  of  our  predictive  models  and  the  validation  data  we 
compiled  from  other  work,  to  prioritize  what  additional  experimental  data  is  necessary  to 
validate  our  models.  We  made  some  of  these  measurements  as  part  of  this  project,  and 
mined  experimental  data  from  other  recent,  SERDP-funded  work  on  the  environmental 
fate  of  insensitive  munitions  compounds. 
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RESULTS  AND  DISCUSSION 


The  main  results  of  this  project  are  presented  in  three  sections  that  correspond  to  the  three 
major  objectives  of  the  proposal.  Each  of  these  sections  is  derived  from  one  or  more  peer- 
reviewed  publications,  which  are  specified  at  the  beginning  of  each  section. 

1.  Objective  1 — Hydrolysis 

The  major  results  of  Objective  1  were  published  in:  Salter-Blanc,  A.  J.,  E.  J.  Bylaska,  J.  J. 
Ritchie,  and  P.  G.  Tratnyek  2013.  Mechanisms  and  kinetics  of  alkaline  hydrolysis  of 
the  energetic  nitroaromatic  compounds  2,4,6-trinitrotoluene  (TNT)  and  2,4- 
Dinitroanisole  (DNAN).  Environ.  Sci.  Technol.  47(13):  6790-6798. 

[10.1 02  l/es30446  It] 

Abstract:  The  environmental  impacts  of  energetic  compounds  can  be  minimized  through 
the  design  and  selection  of  new  energetic  materials  with  favorable  fate  properties. 
Building  predictive  models  to  infonn  this  process,  however,  is  difficult  because  of 
uncertainties  and  complexities  in  some  major  fate-determining  transfonnation  reactions 
such  as  the  alkaline  hydrolysis  of  energetic  nitroaromatic  compounds  (NACs).  Prior  work 
on  the  mechanisms  of  the  reaction  between  NACs  and  OFT  has  yielded  inconsistent 
results.  In  this  study,  the  alkaline  hydrolysis  of  2,4,6-trinitrotoluene  (TNT)  and  2,4- 
dinitroanisole  (DNAN)  was  investigated  with  coordinated  experimental  kinetic 
measurements  and  molecular  modeling  calculations.  For  TNT,  the  results  suggest 
reversible  fonnation  of  an  initial  product,  which  is  likely  either  a  Meisenheimer  complex 
or  a  TNT  anion  formed  by  abstraction  of  a  methyl  proton  by  OFT.  For  DNAN,  the  results 
suggest  that  a  Meisenheimer  complex  is  an  intermediate  in  the  formation  of  2,4- 
dinitrophenolate.  Despite  these  advances,  the  remaining  uncertainties  in  the  mechanisms 
of  these  reactions — and  potential  variability  between  the  hydrolysis  mechanisms  for 
different  NACs — mean  that  it  is  not  yet  possible  to  generalize  the  results  into  predictive 
models  (e.g.,  quantitative  structure-activity  relationships,  QSARs)  for  hydrolysis  of  other 
NACs. 

1.1  Introduction 

Many  energetic  compounds  are  nitroaromatic  compounds  (NACs).  The  most  familiar 
member  of  this  group  is  2,4,6-trinitrotoluene  (TNT),  which  has  long  been  used  in  a 
variety  of  energetic  materials.  Also  included  is  a  shock-insensitive  alternative  to  TNT, 
2,4-dinitroanisole  (DNAN)  [7],  which  is  used  in  modem  munitions  formulations  such  as 
PAX-21  [5],  While  TNT  and  DNAN  comprise  a  subfamily  of  substituted  nitrobenzenes, 
other  NACs  used  in  emerging  munitions  formulations  include  substituted  heterocycles, 
such  as  5-nitro-l,2,4-triazol-3-one  (NTO) — another  insensitive  replacement  for  TNT  that 
is  used  in  emerging  munitions  fonnulations  [7] — and  polyaromatic  compounds.  A 
representative  sample  of  energetic  NACs  is  shown  in  Figure  2;  these  and  many  more 
have  been  described  previously  [ 1 ,  9]. 
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Figure  2.  A  selection  of  energetic  NACs  included  or  considered  for  use  in  munitions 
formulations. 


The  main  abiotic  transformation  processes  influencing  the  environmental  fate  of 
NACs  in  groundwater  are  reduction  and  alkaline  hydrolysis  \10\.  While  the  pathways  and 
kinetics  associated  with  the  reduction  of  NACs  have  been  studied  extensively  (e.g.,  [11- 
14]),  the  pathways  associated  with  alkaline  hydrolysis  are  not  well  characterized,  despite 
considerable  proof-of-concept  testing  of  this  process  for  engineered  remediation  of  TNT- 
contaminated  wastewaters  [15-23].  In  the  case  of  TNT,  there  are  significant 
inconsistencies  in  the  experimental  data,  which  likely  arise  from  the  sensitivity  of  the 
reaction  to  the  nature  of  the  solvent  [24,  25],  difficulty  observing  products  in  water  due  to 
their  poor  solubility  [19],  and  the  influence  of  the  concentration  ratio  between  TNT  and 
base  on  product  distribution  [24,  25],  difficulty  observing  products  in  water  due  to  their 
poor  solubility  [19].  Previous  computational  studies  of  TNT  reaction  energetics  have  also 
failed  to  definitively  detennine  the  pathways  of  its  reaction  with  OH-  [26-28]. 

This  incomplete  understanding  of  the  mechanisms  of  TNT  hydrolysis  in  natural  and 
engineered  systems  inhibits  the  prediction  of  degradation  properties  of  future  energetic 
NACs  through  the  use  of  correlations  calibrated  with  data  detennined  fully  in  silico.  Our 
goal  under  this  objective  was  to  clarify  the  mechanisms  of  TNT  degradation  by  OH-  in 
water  by  taking  a  combined  approach  that  emphasizes  reconciliation  of  both  experimental 
and  computational  data.  An  analysis  of  the  kinetics  of  TNT  disappearance  in  the  presence 
of  OH“  is  presented  and  used  to  detennine  the  experimental  activation  free  energies 
(AG:)  and  reaction  free  energies  (A  G,Xn)  for  the  processes.  Also  presented  arc  AG?  and 
AGrxn  values  for  possible  mechanisms  of  TNT  degradation  detennined  using  molecular 
modeling.  The  measured  and  modeled  results  are  then  compared  to  assess  which 
mechanisms  predominate.  Also  reported  is  an  analogous  analysis  for  DNAN,  which  is 
structurally  similar  to  TNT. 

1.2  Methods 

Method  details  are  more  fully  documented  in  Salter-Blanc  et  al.  [29]. 

Batch  experiments.  Batch  experiments  were  canied  out  in  20-mL  amber  VOA  vials 
capped  with  Teflon-lined  silicon  septa  (SUN  Sri;  Rockwood,  TN).  Each  vial  initially 
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contained  20  mL  phosphate  buffer  (50  mM)  at  pH  1 1.0,  1 1.7,  or  12.0.  Vials  were 
temperature  equilibrated  in  water  baths  or  in  a  cold  room  to  1 1.0,  25.0,  40.0,  55.0,  or  65.0 
°C.  After  temperature  equilibration,  TNT  or  DNAN  was  introduced  to  the  reaction  vial 
by  injecting  200  pL  of  a  10-g/L  stock  solution  prepared  in  acetonitrile  (reaction 
initiation).  Following  introduction  of  TNT  or  DNAN,  the  vials  were  shaken  by  hand  for 
~1  min  (in  the  water  bath)  to  insure  proper  mixing.  1-mL  aliquots  were  removed  at 
specified  times  and  quenched  by  mixing  with  an  equal  volume  of  acidified  acetonitrile. 
Quenched  aliquots  were  analyzed  by  high-pressure  liquid  chromatography  (HPLC). 

Molecular  Modeling.  In  this  study,  the  solution  phase  AGrxn  and  AG*  were  directly 
calculated  from  gas-phase  reaction  energy,  entropy,  electronic  structure  calculations, 
continuum  solvation  models,  and  gas-phase  entropy  estimates.  All  calculations  were 
perfonned  using  the  NWChem  program  suite  [50].  The  electronic  structure  calculations 
were  either  perfonned  using  density  functional  theory  (DFT)  [57]  with  the  B3LYP  [32, 
55]  exchange  correlation  potential,  or  with  second-order  Moller-Plesset  perturbation 
theory  (MP2)  [54].  In  all  cases,  the  6-3 1 1++G(2d,2p)  basis  set  was  used  (obtained  from 
the  Extensible  Computational  Chemistry  Environmental  Basis  Set  Database  [55]).  The 
solvation  energies  were  estimated  using  the  self-consistent  reaction  field  theory  of  Klamt 
and  Schiiurmann  (COSMO)  [36],  with  the  cavity  defined  by  a  set  of  overlapping  atomic 
spheres  with  radii  suggested  by  Stefanovich  and  Truong  [37].  The  dielectric  constant  of 
water  used  for  all  of  the  solvation  calculations  was  78.4.  The  solvent  cavity  discretization 
was  generated  from  the  surfaces  of  non-overlapping  spheres  that  were  discretized  by  an 
iterative  refinement  of  triangles  starting  from  a  regular  octahedron. 

The  geometries  and  harmonic  frequencies  for  the  reactants,  products  and  transition 
states  were  consistently  optimized  using  the  B3LYP  calculation  with  COSMO. 
Transitions  states  were  detennined  for  the  Mei sen hei trier  and  proton  abstraction  reactions 
by  perfonning  constrained  optimizations  at  a  series  of  C-OH  and  CH-OH  bond  distances. 
The  virtual  entropies  for  each  compound  were  estimated  using  formulas  derived  from 
statistical  mechanics  that  are  broken  into  translational,  rotational,  and  frequency  terms.  In 
addition  to  the  COSMO  correction,  cavitation,  and  dispersion  contributions  to  the 
solvation  energy  were  added  a  posteriori  using  empirically-derived  expressions  that 
depend  only  on  the  solvent  accessible  surface  area.  In  this  study,  we  used  the 
parameterized  formula  given  by  Sitkoff  et  al.  [55]. 

1.3  Results 

1 .3.1  Formulation  of  hydrolysis  reaction  pathways/mechanisms  and  mining 
previously-published  data  on  the  reactants  and  products  (Task  1.1) 

1. 3.1.1  Mechanism  Determination 

Cyclic  and  cage-cyclic  nitramines.  The  most  highly-studied  cyclic  nitramine,  with 
respect  to  hydrolysis  reaction  mechanisms,  is  RDX.  Aqueous  alkaline  hydrolysis  of  RDX 
occurs  through  a  concerted  two-step  process  involving  proton  abstraction  by  hydroxide 
and  loss  of  a  nitrite  ion  (NCfi-)  [39]  (E2  elimination  reaction).  Initial  attack  by  hydroxide 
is  rate-limiting  [40,  41]  and  the  reaction  is  second  order  with  respect  to  hydroxide 
concentration  [39,  41-44].  Denitration  is  followed  by  ring  cleavage  and  spontaneous 
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decomposition  to  N2O,  HCHO,  and  4-NDAB  [ 45] .  Identical  products  are  observed 
following  hydrolysis  of  HMX  and  MNX,  which — like  RDX — are  observed  to  produce 
one  equivalent  of  NCk”  in  the  reaction  process  [45].  This  supports  earlier  suggestions  that 
cyclic  nitramines  other  than  RDX  also  react  via  an  E2  elimination  reaction  [46].  One 
theoretical  study  suggests  that  an  SnI  hydroxide  ion  substitution  of  RDX  is  more 
energetically  favorable  than  the  E2  elimination  reaction  [47],  but  this  suggestion  is  in 
disagreement  with  experimentally  observed  products  [45]. 

Structurally  similar  to  the  cyclic  nitramines  is  CL-20,  a  cage-cyclic  (or  polycyclic) 
nitramine.  Kinetic  studies  of  CL-20  show  second-order  behavior  with  respect  to  CL-20 
and  hydroxide  concentrations  [45,  48,  49].  One  study  suggests  the  reaction  proceeds  by 
initial  denitration — similarly  to  the  monocyclic  nitramines — however,  the  reaction 
produces  two  equivalents  of  NO2”  and  different  final  products  than  were  observed  in  the 
monocyclic  experiments.  It  has  also  been  suggested  that  this  reaction  proceeds  differently 
under  high  or  low  concentrations  of  hydroxide,  producing  products  that  have  different 
susceptibilities  to  photodegradation  [47].  Comparison  of  experimental  FTIR  spectra  for 
CL-20  hydrolysis  and  DFT  theorized  spectra  for  1,5-  and  1,7-  dihy  dr  odiimidazo  [  7, 5-h :  7 ' 
,5'-<?]pyrazine  suggest  CL-20  hydrolysis  produces  a  mixture  of  these  tautomers  through 
cleavage  of  the  C-C  bond  attaching  the  two  5-membered  rings  [50]. 

Nitroaromatic  compounds.  Most  relevant  environmental  literature  on  the  alkaline 
hydrolysis  of  nitroaromatic  compounds  focuses  on  TNT.  The  reaction  of  TNT  with 
hydroxide  is  second  order  with  respect  to  TNT  and  hydroxide  concentrations  [16,  19]. 
The  pathways  of  TNT  hydrolysis  are  complex — intennediate  and  product  formation 
appear  to  be  affected  by  pH  [23],  the  TNT/hydroxide  concentration  ratio  [27],  and  the 
reaction  media  [19].  Prior  studies  led  to  three  main  hypotheses  regarding  the  first  step  of 
TNT  alkaline  hydrolysis;  these  are:  (/)  formation  of  a  Meisenheimer  complex  [17,  19, 

51],  (ii)  direct  nitro  substitution  [20],  and  (iii)  abstraction  of  a  methyl  proton  [27,  25], 
More  details  are  given  in  the  publication  resulting  from  this  objective  [29].  Alternative 
initial  steps  in  the  reaction  between  TNT  and  OH“  are  shown  in  Figure  3A.  A  limited 
number  of  studies  have  addressed  the  alkaline  hydrolysis  of  2,4-dinitroanisole  (DNAN). 
The  reaction  is  suggested  to  occur  by  substitution  of  the  methoxy  group  to  form  2,4- 
nitrophenol,  which  in  turn  may  react  with  OH“  to  from  2,4-dinitrophenolate  (Figure  3B). 
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Meisenheimer  complex  formation 


NOj 

Proton  abstraction 


Figure  3.  (A)  Alternative  initial  steps  in  the  mechanism  of  reaction  between  TNT  and 
OH“.  (B)  Alternative  initial  steps  in  the  mechanism  of  reaction  between  DNAN  and  OH“. 
Meisenheimer  complex  formation  at  Cl,  which  has  been  hypothesized  based  on 
experimental  observations  [52],  is  shown  in  bold. 


1.3. 1.2  Determination  of  “Best  Values”  of  kinetic  data  from  Literature 

Rate  constants  for  the  base-catalyzed  hydrolysis  of  various  nitramines  were  collected 
from  the  literature.  Reported  values  for  ki,  the  second-order  rate  constant,  were  compared 
in  Arrhenius  plots.  In  some  cases,  ki  was  not  reported  directly,  but  was  calculated  from 
reported  values  of  k0 bs,  the  pseudo-first-order  rate  constant,  measured  at  various  pH’s. 

This  set  of  data  was  used  to  detennine  the  temperature  dependence  of  the  reaction  as 
defined  by  the  Arrhenius  equation  (Equation  1),  where  A  is  the  pre-exponential  factor  and 
Ea  is  the  activation  energy. 


k2=Je 


EJRT 


The  natural  logarithm  of  Equation  1  is  shown  in  Equation  2: 


(i) 


In {k  )  =  — — + ln(^) 

2  R  T  (2) 

Equation  2  was  used  to  fit  ki  vs.  1/T  data  for  RDX  (Figure  4 A),  HMX  (Figure  5A)  and 
CL-20  (Figure  6A)  to  obtain  the  Arrhenius  parameters  A  and  Ea.  In  all  cases,  the  fit 
shows  excellent  agreement  between  the  ki  values  obtained  from  different  sources.  The 
only  data  point  that  appears  to  deviate  is  the  value  for  RDX  calculated  from  data  reported 
by  Hwang,  et  al.  [ 43\  (Figure  4A).  This  data  point  was  excluded  from  the  final  fit.  Also 
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shown  is  a  plot  for  TNT  (Figure  7A).  Agreement  between  the  values  in  this  figure  is  not 
as  good  as  seen  with  the  other  compounds. 

The  Arrhenius  parameters  obtained  from  fitting  the  data  in  Figure  4A,  Figure  5A, 
and  Figure  6A  were  used  to  predict  the  dependence  of  the  observed  rate  constants  (k0 bs) 
on  pH  and  temperature  by  using  the  relationship  k0 bs  =  ki  [OH~]  and  Equation  5.  The 
results  are  shown  in  Figure  4B,  Figure  5B,  and  Figure  6B  with  diagonal  lines.  Also  on 
these  plots  are  the  k0 bs  data  used  by  Hoffstetter,  et  al.  [39],  Heilmann,  et  al.  [42],  Hwang, 
et  al.  [43],  Monteil-Rivera,  et  al.  [44],  Karakaya,  et  al.  [48],  and  Santiago,  et  al.  [49]  to 
obtain  the  ki  values  shown  in  Figure  4A,  Figure  5A,  and  Figure  6A.  Additional  k0 bs  data 
from  Balakrishnan,  et  al.  [45]  are  also  shown,  although  these  data  were  collected  under 
heterogeneous  conditions,  whereas  the  other  data  were  collected  under  homogeneous 
conditions.  With  the  exception  of  the  Balakrishnan,  et  al.  [45]  data,  the  figures  show 
excellent  agreement  between  the  values  reported  by  the  different  sources  and  the 
calculated  Arrhenius  parameters.  This  agreement  supports  the  use  of  these  parameters  to 
calculate  “best  values”  for  the  hydrolysis  rate  constants  to  be  used  in  correlation  analysis. 


Temperature  (X) 

80  60  40  20 


1/T  (1/K) 


pH 


Figure  4.  RDX.  (A)  Arrhenius  plot  displaying  reported  and  calculated  data  from 
designated  sources  [39,  42-44],  Data  was  fit  to  the  displayed  linear  model  to  obtain  the 
collective  Arrhenius  parameters.  The  data  point  from  Hwang,  et  al.  [43]  was  not  included 
in  the  regression  analysis.  (B)  Dependence  of  kD bS  on  temperature  and  pH.  Isothermal 
trends  for  differing  temperatures  calculated  using  the  Arrhenius  parameters  determined 
in  part  A  of  this  figure  are  shown  with  diagonal  lines.  Reported  ka bS  values  from 
designated  sources  are  also  displayed  [39,  42-45], 
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Temperature  (°C) 
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Figure  5.  HMX.  (A)  Arrhenius  plot  displaying  reported  ta  data  from  designated  sources 
[42,  44],  Data  was  fit  to  the  displayed  linear  model  to  obtain  the  collective  Arrhenius 
parameters.  (B)  Dependence  of  tabs  on  temperature  and  pH.  Isothermal  trends  for 
differing  temperatures  calculated  using  the  Arrhenius  parameters  determined  in  part  A  of 
this  figure  are  shown  with  diagonal  lines.  Reported  tabs  values  from  designated  sources 
are  also  displayed  [42,  44]. 
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Figure  6.  CL-20.  (A)  Arrhenius  plot  displaying  reported  and  calculated  ta  data  from 
designated  sources  [48,  49],  Data  was  fit  to  the  displayed  linear  model  to  obtain  the 
collective  Arrhenius  parameters.  (B)  Dependence  of  tabs  on  temperature  and  pH. 
Isothermal  trends  for  differing  temperatures  calculated  using  the  Arrhenius  parameters 
determined  in  part  A  of  this  figure  are  shown  with  diagonal  lines.  Reported  tabs  values 
from  designated  sources  are  also  displayed  [45,  48,  49]. 
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Figure  7.  TNT.  Arrhenius  plot  displaying  reported  and  calculated  fe  data  from 
designated  sources  [16,  19}.  Data  was  fit  to  the  displayed  linear  model  to  obtain  the 
collective  Arrhenius  parameters. 

1.3. 1.3  Organization  of  Data  in  a  Chemically-lntelligent  Database 

In  addition  to  mining  mechanism  and  rate  constant  data  from  the  literature,  another  aspect 
of  Task  1.1  was  the  development  of  a  chemoinformatic  database  for  storing  this  and  other 
data.  After  inspecting  chemically  intelligent  database  packages,  we  decided  to  use  Instant 
JChem  by  ChemAxon  (Budapest,  Hungary).  This  software  allows  for  storage  and  query 
of  numerical,  textual,  and  structural  data  and  for  the  development  of  relational  databases. 
This  will  be  very  useful  for  storing  complex  data  and  meta  data  and  for  sorting  by 
mechanistically-relevant  structural  moieties.  The  ChemAxon  tools  also  include  certain 
property  predictors  such  as  a  pKa  calculator,  which  may  be  used  in  determining 
descriptors  for  QSAR  development.  Screen  shots  of  the  database  are  shown  in  Figure  8 
to  Figure  10. 

The  decision  to  use  Instant  JChem  was  partly  motivated  by  the  adoption  of 
ChemAxon  tools  by  Dr.  Eric  Weber’s  group  of  the  EPA,  Athens  (for  his  SERDP  project, 
but  also  for  broader  goals  of  EPA’s  “CompTox”  programs.  Dr.  Weber’s  group  has  been 
working  on  developing  an  environmental  fate  simulator  for  novel  and  proposed 
munitions  compounds.  The  co-adoption  of  ChemAxon  tools  by  Dr.  Weber’s  group  and 
our  group  facilitated  collaboration  and  the  integration  of  our  results  (predicted  using 
QSARs  calibrated  with  fully  in  silico  methods)  into  the  proposed  fate  simulator. 
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Figure  8.  Screen  shot  of  a  spreadsheet-like  view  of  the  database,  showing  basic 
information  about  each  compound. 
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Figure  9.  “Form  view”  page  for  an  individual  compound  (RDX)  with  data  linked  from  the 
database  component  shown  in  Figure  8  and  from  another  component  containing  kinetic 
data  mined  from  the  literature. 
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Figure  10.  Results  of  a  query  for  a  structural  moiety  involved  in  the  mechanism 
described  above  for  the  hydrolysis  of  RDX  and  related  compounds.  It  is  likely 
compounds  in  this  set  hydrolyze  by  a  similar  mechanism,  so  it  will  be  logical  to  group 
these  compounds  in  the  same  QSAR. 

1 .3.2  Calculation  of  target/response  variable  data  for  training  set  compounds 
with  high-level  theory  (Task  1 .2). 

As  reviewed  for  Tasks  1.1,  there  are  several  hydrolysis  mechanisms  for  explosive 
compounds.  The  three  most  likely  are  substitution, 

R-NCh  +  OH-  —  R-OH  +  N02"  (3) 

Meisenheimer  addition  where  the  hydroxyl  attaches  to  a  neighboring  carbon  atom  in  an 
aromatic  ring, 

Ar-N02  +  OH"  —  (0H)-Ar-N02"  (4) 

and  the  highly  favorable  [1-elimination  reaction  for  nitramine  explosives. 

R-CH-N-N02  +  OH-  —  R-C=N  +  H20  +  N02"  (5) 

Examples  of  our  calculated  reaction  free  energies  for  a  few  of  the  nitroaromatics  are 
shown  in  Table  1  and  Table  2.  In  general  for  the  nitroaromatic  compounds  it  was  found 
that  the  overall  thermodynamics  of  the  substitution  reactions  was  -20  kcal/mol  more 
favorable  than  the  Meisenheimer  reactions  for  the  same  compound.  It  was  also  found 
that  Meisenheimer  complexes  were  only  thermodynamically  stable  when  the 
nitroaromatic  compounds  contained  3  nitro  groups  (e.g.  TNT,  TNB,  tetryl). 
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Table  1.  Free  energies  of  reaction  in  kcal/mol  for  substitution  reactions  involving 
nitroaromatic  compounds  calculated  at  PBE,  B3LYP,  and  PBEO  theories  with  COSMO 
solvation  model. 


Substitution  Reactions 

PBE 

B3LYP 

PBEO 

TNT  +  OH"  ->  TNT-2-OH  +  N02“ 

-30.9 

-33.5 

-33.3 

TNT  +  OH"  ->  TNT-4-OH  +  N02“ 

-28.4 

-30.4 

-30.9 

TNT  +  OH"  ->  TNT-1-OH-2-CH3  +  N02" 

-28.8 

-32.1 

-31.1 

TNT  +  OH"  ->  TNT-3-OH-4-H  +  N02" 

-28.7 

-30.0 

-30.1 

TNT  +  OH"  ->  TNT-3-OH-2-H  +  N02“ 

-29.5 

-34.1 

-33.7 

NB  +  OH"  ->  NB-1-OH  +  N02“ 

-21.4 

-23.7 

-23.8 

TNB  +  OH"  -h.  TNB-1 -OH-2, 3-N02  +  N02“ 

-24.6 

-25.7 

-25.1 

TNB  +  OH"  -h.  TNB-1 -OH-2, 4-N02  +  N02“ 

-34.6 

-36.7 

-35.9 

TNB  +  OH"  ->■  TNB-1 -OH-2, 5-N02  +  N02" 

-31.4 

-33.0 

-32.5 

TNB  +  OH"  -h.  TNB-1 -OH-2, 6-N02  +  N02" 

-27.4 

-29.1 

-28.5 

TNB  +  OH"  -h.  TNB-1 -OH-3, 4-N02  +  N02" 

-27.8 

-30.0 

-30.1 

TNB  +  OH"  -h.  TNB-1 -OH-3, 5-N02  +  N02“ 

-31.1 

-34.5 

-33.6 

DNA  +  OH"  -h.  DNA-2-OH"  +  N02" 

-25.1 

-28.2 

-28.1 

TATB  +  OH"  -»  TATB-2-OH  +  N02" 

-24.8 

-28.9 

-28.1 

Tetryl  +  OH-  — tetryl-2-OH  +  N02“ 

-33.2 

-36.3 

-36.4 

Tetryl  +  OH"  — » tetryl-4-OH  +  N02" 

-34.2 

-35.4 

-34.8 

DNAN  +  OH"  -h.  DNAN-2-OH  +  N02“ 

-24.4 

-27.8 

-26.6 

DNAN  +  OH"  -h.  DNAN-4-OH  +  N02“ 

-14.1 

-17.6 

-16.0 

Table  2.  Free  energies  of  reaction  in  kcal/mol  for  Meisenheimer  addition  reactions 
involving  nitroaromatic  compounds  calculated  at  PBE,  B3LYP,  and  PBEO  theories  with 
COSMO  solvation  model. 


Meisenheimer  Reactions 

PBE 

B3LYP 

PBEO 

TNT  +  OH"  — >  TNT-1-OH" 

-9.6 

-3.8 

-10.1 

TNT  +  OH"  ->  TNT-3-OH" 

-7.9 

-1.7 

-8.5 

NB  +  OH"  ->  NB-2-OH" 

13.7 

19.7 

13.5 

NB  +  OH"  ->  NB-3-OH" 

14.2 

14.9 

14.3 

NB  +  OH"  -»  NB-4-OH" 

15.0 

14.2 

14.0 

TNB  +  OH-  — >  TNB-2-OH- 

-15.8 

-11.2 

-16.4 

DNA  +  OH-  — >  DNA-1-OH- 

-0.5 

5.1 

-0.9 

DNA  +  OH"  -h.  DNA-3-OH" 

2.9 

9.1 

3.9 

Tetryl  +  OH"  — ►  tetryl-1-OH" 

-17.1 

-10.7 

-19.2 

Tetryl  +  OH"  — >  tetryl-3-OH" 

-9.0 

-7.8 

-3.6 

Tetryl  +  OH"  — >  tetryl-4-OH" 

-15.9 

-12.7 

-14.2 

DNAN  +  OH"  DNAN-1-OH- 

0.7 

5.3 

0.9 

DNAN  +  OH"  -h.  DNAN-3-OH" 

19.1 

16.6 

DNAN  +  OH"  -h.  DNAN-5-OH" 

19.0 

16.2 

19.9 

In  our  initial  calculations,  we  attempted  to  determine  the  activation  barriers  of  the 
substitution  reactions  of  nitroaromatics  by  OH"  using  a  two-step  process,  where  we  first 
calculated  the  transition- state  and  activation  barrier  for  the  hydrolysis  reaction  in  the  gas- 
phase.  Followed  by  calculation  of  solvent  effects  the  activation  barriers  using  continuum 
solvation  models  and  QM/MM  models.  However,  it  turns  out  that  the  hydrolysis  reaction 
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does  not  have  a  activation  barrier  in  the  gas-phase,  because  its  overall  reaction  energy  is 
extremely  favorable  (e.g.  CfitUNCFUg)  +  OH“(g)  — >  C6H50H(g)  +  NC>2-(g),  Erxn(B3LYP)  = 
-50  kcal/mol). 

Given  difficulties  with  this  approach,  we  had  to  use  computationally  more  expensive 
strategies  to  determine  these  activation  barriers.  The  approach  that  worked  best  was  to 
map  out  the  two-dimensional  reaction  energy  surface  in  solution  (i.e.  with  a  continuum 
model)  as  a  function  of  the  R-NO2  and  R-OH  distances  (Figure  11).  This  approach  was 
considerably  more  expensive  then  our  initial  approach  as  it  requires  ~500  constrained 
geometry  optimizations  to  compute  each  reaction  energy  surface.  For  example  at  the  DFT 
level,  it  takes  -150K  CPU  hours  (0.73  CPU  years)  to  calculate  a  surface  for  the  reaction 
of  TNT  +  OH-.  Using  this  approach,  we  have  detennined  the  reaction  activation  barriers 
for  TNT,  NB,  1,4-DNB,  and  1,3-DNB  (Table  4). 


Figure  11.  (a)  Two  dimensional  potential  energy  surface  for  the  reaction,  C6H5N02(aq)  + 
OH“(aq)  -*•  C6H5OH(aq)  +  NO2-  (aq),  calculated  at  the  B3LYP/6-31 1++G(3p,3d)  //COSMO 
level.  The  reaction  barrier  for  this  potential  energy  surface  is  Eact=  +31  kcal/mol.  (b) 
Illustration  of  an  AIMD/MM  simulation  for  the  C6H5N02(aq)  +  OH~(aq)  — ►  C6H5OH(aq)  +  NO2- 
(aq>  reaction  and  the  two  dimensional  free  energy  surface  of  the  reactant  well  obtained 
from  the  AIMD/MM  metadynamics  simulation.  The  reaction  barrier  for  this  potential 
energy  surface  is  Eact=  +30  kcal/mol. 


Table  4.  Activation  barriers  in  kcal/mol  for  TNT,  NB,  1,4-DNB,  and  1,3-DNB  from  two 
dimensional  potential  energy  surfaces  calculated  at  the  B3LYP/6- 
31 1  ++G(3p,3d)//COSMO  level. 


Compound 

Eact 

TNT 

20  kcal/mol 

NB 

31  kcal/mol 

1,4-DNB 

15  kcal/mol 

1,3-DNB 

23  kcal/mol 
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We  have  also  been  using  AIMD/MM  (33)  free  energy  approaches  to  calculate 
potential  energy  surfaces  for  OH“  substitution  reactions.  In  these  types  of  simulations,  it 
is  sufficient  to  compute  AG  and  AG+  of  reactions  as  just  the  Helmholtz  free  energy,  A A, 
and  activation  barrier,  A/4+,  because  the  additional  terms  PAV and  PAV*  in  AG  and  AG+ 
respectivey  are  usually  very  small.  We  used  the  metadynamics  free  energy  sampling 
method  with  an  AIMD/MM  model  (Figure  11)  to  detennine  the  A A  and  A A*  for  the 
various  OHYNO2  exchange  processes  in  nitro  containing  compounds.  Metadynamics  is 
a  powerful,  non-equilibrium  molecular  dynamics  method  that  accelerates  the  sampling  of 
the  multidimensional  free  energy  surfaces  of  chemical  reactions.  For  the  OH“  substitution 
the  reactions,  these  simulations  can  be  perfonned  on  a  modest  workstation  in  about  a 
week. 

1 .3.3  Calculate  descriptor  variable  data  for  training  set  compounds  using  lower 
level  theory  (Task  1 .3) 

Energy-related  descriptor  variables  such  as  Er,  Ehomo,  Elumo,  and  Egap  were  calculated 
using  relatively-accessible  methods.  Calculations  were  perfonned  on  a  desktop  computer 
(MacMini,  2.66  GHz  Intel  Core  2  Dual  processer)  using  Gaussian09  with  the  GausView5 
interface.  These  calculations  were  performed  using  density  functional  theory  (DFT)  with 
the  B3LYP  functional  and  6-31G(d)  basis  set.  These  calculations  were  successful,  but  the 
efficient  optimization  of  larger  structures  (e.g.,  CL-20)  required  greater  computational 
resources  than  provided  by  the  desktop  system.  In  order  to  provide  a  less  computationally 
expensive  alternative,  a  recommended  “research-level”  model  chemistry  was  also 
employed.  In  this  method,  B3LYP/6-31G  energy  calculations  were  perfonned  on 
structures  optimized  using  Hartree-Fock  (HF)  theory  (HF/3-21G(d)). 

A  comparison  of  Elumo  values  calculated  from  DFT  and  HF  optimized  structures  is 
shown  in  Figure  12.  The  overall  correlation  (regression  line)  shows  a  near  1:1 
relationship,  with  a  slope  of  0.996,  intercept  of  0.024,  and  scatter  about  the  regression 
line.  The  significance  of  these  variations  (especially  scatter)  with  respect  to  QSAR 
formulation  and  is  discussed  below. 
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Figure  12.  Comparison  of  DFT  calculations  of  Elumo  from  structures  optimized  using  HF 
theory  and  DFT. 

1 .3.4  Correlation  Analysis  and  fitting  QSARs  (Task  1 .4) 

Two  groups  of  compounds  were  selected  for  QSAR  prediction  of  hydrolysis  rate 
constants  based  on  structural  similarities  likely  to  lead  to  similar  mechanisms.  These 
groupings  are  “nitramine-like”  compounds  (i.e.,  with  a  C-N(-NC>2)  moiety)  and 
nitroaromatic  compounds,  with  our  main  focus  being  on  the  nitro  aromatic  compounds. 
The  reaction  of  these  compounds  with  hydroxide  in  water  is  well  understood  [57]  and 
there  are  a  variety  of  existing  and  emerging  munitions  compounds  that  fit  this 
description,  as  well  as  a  number  of  model  compounds. 

Although  we  had  a  large  number  of  nitroaromatic  compounds  in  our  database  (>60), 
we  focused  on  compounds  with  a  nitrobenzene  backbone.  We  designated  a  handful  of 
these  compounds  as  “top-priority”.  The  selection  of  these  compounds  was  based  on 
factors  including  novelty,  designation  as  a  SEDP  priority  compound,  structural  similarity 
to  novel  or  priority  compounds,  current  usage,  previous  publication  history,  or 
availability  for  use  in  bench-top  experiments  (munitions  and  model  compounds).  We  also 
further  organized  the  list  into  “Priority  I”  and  “Priority  II”  compounds,  where  the  Priority 
I  group  contains  munitions  compounds  and  important/available  model  compounds,  and 
the  Priority  II  group  contains  a  number  substituted  nitrobenzenes  that  have  been  well 
studied  by  the  environmental  chemistry  community  [11,  12,  53-55]. 

Our  current  list  is  shown  in  Table  3.  We  also  plan  to  include  2,6-diamino-3,5- 
dinitropyrazine- 1  -oxide  (LLM-105)  in  our  study,  although  it  is  not  a  substituted 
nitrobenzene,  because  of  its  importance  as  an  emerging  munitions  compound. 
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Table  3.  List  of  priority  compounds  for  QSAR  development. 


Priority  1 

Priority  II 

2,4,6-Trinitroluene  (TNT) 

1,4-Dinitrobenzene  (1-4-DNB) 

2,4-Dinitroanisole  (DNAN) 

2,4-Diamino-6-nitrotoluene  (2,4-DANT) 

Nitrobenzene  (NB) 

2,6-Diamino-4-nitrotoluene  (2,6-DANT) 

1,3-Dinitrobenzene  (1,3-DNB) 

2-Methylnitrobenzene  (2-CH3-NB) 

2-Aminodinitrobenzene  (2-ADNT) 

2-Acetylnitrobenzene  (2-COCH3-NB) 

3,5-Dinitroanaline  (DNA) 

2-Nitroanaline  (2-NH2-NB) 

1, 3, 5-Triamino-2, 4, 6-trinitrobenzene  (TATB) 

3-Methylnitrobenzene  (3-CH3-NB) 

Trinitrobenzene  (TNB) 

3-Acetylnitrobenzene  (3-COCH3-NB) 

Tetryl 

3-Nltroanaline  (3-NH2-NB) 

Musk  ketone 

4-COCH3-NB 

Musk  xylene 

4-Amino-2,6-dinitrotoluene  (4-ADNT) 

1 .3.5  Verification  of  QSARs  with  experimental  values  (Task  1 .5) 

We  measured  hydrolysis  rates  for  verifying  the  predictability  of  potential  QSARs.  As  a 
first  step,  we  determined  rate  constants  for  the  hydrolysis  of  TNT.  This  has  allowed  us  to 
compare  our  values  to  those  reported  in  the  literature,  and  thereby  check  of  the  accuracy 
of  our  methods.  We  will  also  detennine  rate  constants  for  the  hydrolysis  of  other 
nitroaromatic  compounds  with  a  nitrobenzene  backbone  (both  energetic  materials  and 
model  compounds).  The  results  for  the  TNT  experiments  are  shown  in  Figure  13  along 
with  data  reported  in  the  literature  [16,  19].  There  is  some  scatter  in  the  data.  Part  of  the 
scatter  in  the  data  is  likely  due  to  the  fact  that  Mills  et  al.  [19]  fit  their  concentration  vs. 
time  data  to  a  reversible  first-order  model  (which  is  appropriate  for  this  reaction), 
whereas  Emmrich  [16]  neglected  to  take  into  account  the  reversibility  of  the  reaction.  We 
have  taken  into  account  the  reversibility  of  the  reaction,  yet  our  data  seems  to  fall  in 
between  the  data  gathered  by  Mills  et  al.  and  Emmrich.  The  reason  for  this  is  not  clear  at 
this  time. 
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Figure  13.  (A)  Arrhenius  plot  and  (B)  Eyring  plot  showing  second  order  rate  constants 
for  alkaline  hydrolysis  of  TNT.  Both  new  data  and  data  reported  in  the  literature  are 
shown.  Arrhenius  parameters  as  well  as  the  enthalpy  and  entropy  of  the  transition  state 
(4/-/*  and  AS*)  were  obtained  from  fitting  the  data  to  the  appropriate  model.  The  free 
energy  for  the  transition  state  (AG*)  at  25  °C  was  calculated  from  AH*  and  AS*. 

1.4  Conclusions 

Previous  efforts  have  failed  to  unambiguously  define  mechanisms  for  the  reactions  of 
TNT  and  DNAN  with  OH“,  especially  in  the  case  of  TNT  where  product  characterization 
has  been  particularly  challenging.  The  experimental  and  computational  observations 
reported  here  provide  insight  into  these  mechanisms,  although  some  ambiguity  remains, 
especially  in  the  case  of  TNT.  Because  of  this  ambiguity,  it  is  uncertain  whether  TNT  and 
DNAN  react  by  the  same  mechanism.  The  possible  difference  in  mechanisms  means  that 
predicting  the  reaction  mechanism  for  one  based  on  the  other  may  lead  to  unreliable 
predictions  of  environmental  fate.  This,  along  with  uncertainties  in  the  consistency  of  the 
calculated  results  with  experimental  values,  presents  a  challenge  for  developing  QSARs 
calibrated  “ in  silico ”  that  predict  the  hydrolysis  behavior  of  the  diverse  range  of  energetic 
NACs. 

2.  Objectives  2  and  3 — Nitro  Reduction 

The  major  results  of  Objective  2/3  were  published  in  two  manuscripts: 

(1)  Salter-Blanc,  A.  J.,  E.  J.  Bylaska,  H.  Johnston,  and  P.  G.  Tratnyek  2015.  Predicting 
reduction  rates  of  energetic  nitroaromatic  compounds  using  calculated  one-electron 
reduction  potentials.  Environ.  Sci.  Technol.  49(6):  3778-3786.  [10.1021/es505092s] 

Abstract.  The  evaluation  of  new  energetic  nitroaromatic  compounds  (NACs)  for  use  in 
green  munitions  formulations  requires  models  that  can  predict  their  environmental  fate. 
Previously  invoked  linear  free  energy  relationships  (LFER)  relating  the  log  of  the  rate 
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constant  for  this  reaction  (log(k))  and  one-electron  reduction  potentials  for  the  NAC 
(El nac)  nonnalized  to  0.059  V  have  been  re-evaluated  and  compared  to  a  new  analysis 
using  a  (non-linear)  free-energy  relationship  (FER)  based  on  the  Marcus  theory  of  outer- 
sphere  electron  transfer.  For  most  reductants,  the  results  are  inconsistent  with  simple  rate 
limitation  by  an  initial,  outer-sphere  electron  transfer,  suggesting  that  the  strong  linear 
correlation  between  log(k)  and  E'nac:  is  best  regarded  as  an  empirical  model.  This 
correlation  was  used  to  calibrate  a  new  quantitative  structure-activity  relationship 
(QSAR)  using  previously  reported  values  of  log(k)  for  non-energetic  NAC  reduction  by 
Fe(II)  porphyrin  and  newly  reported  values  of  E1NAC  detennined  using  density 
functional  theory  at  the  M06-2X/6-3 1 1++G(2d,2p)  level  with  the  COSMO  solvation 
model.  The  QSAR  was  then  validated  for  energetic  NACs  using  newly  measured  kinetic 
data  for  2,4,6-trinitrotoluene  (TNT),  2,4-dinitrotoluene  (2,4-DNT),  and  2,4-dinitroanisole 
(DNAN).  The  data  show  close  agreement  with  the  QSAR,  supporting  its  applicability  to 
other  energetic  NACs. 

(2)  Bylaska,  E.  J.,  A.  J.  Salter-Blanc,  and  P.  G.  Tratnyek.  2011.  One-electron  reduction 
potentials  from  chemical  structure  theory  calculations.  In:  P.  G.  Tratnyek,  T.  J. 
Grundl,  and  S.  B.  Haderlein  (ed.),  Aquatic  Redox  Chemistry.  ACS  Symposium  Series, 
American  Chemical  Society,  Washington,  DC,  Vol.  1071,  pp  37-64. 

Abstract:  Many  redox  reactions  of  importance  in  aquatic  chemistry  involve  elementary 
steps  that  occur  by  single-electron  transfer  (SET).  This  step  is  often  the  first  and  rate 
limiting  step  in  redox  reactions  of  environmental  contaminants,  so  there  has  been  a  great 
deal  of  interest  in  the  corresponding  one-electron  reduction  potentials  (El).  Although  E  1 
can  be  obtained  by  experimental  methods,  calculation  from  first-principles  chemical 
structure  theory  is  becoming  an  increasingly  attractive  alternative.  Sufficient  data  are  now 
available  to  perfonn  a  critical  assessment  of  these  methods — and  their  results — for  two 
types  of  contaminant  degradation  reactions:  dehalogenation  of  chlorinated  aliphatic 
compounds  (CACs)  and  reduction  of  nitro  aromatic  compounds  (NACs).  Early  datasets 
containing  Ehs  for  dehalogenation  of  CACs  by  dissociative  SET  contained  a  variety  of 
errors  and  inconsistencies,  but  the  preferred  datasets  show  good  agreement  between 
values  calculated  from  thermodynamic  data  and  quantum  mechanical  models.  All  of  the 
datasets  with  Eu  s  for  reduction  of  NACs  by  SET  are  relatively  new,  were  calculated  with 
similar  methods,  and  yet  yield  a  variety  of  systematic  differences.  Further  analysis  of 
these  differences  is  likely  to  yield  computational  methods  for  E1  ’s  of  NAC  nitro 
reduction  that  are  similar  in  reliability  to  those  for  CAC  dechlorination.  However, 
comparison  of  the  E 1  data  compiled  here  with  those  calculated  with  a  more  universal 
predictive  model  (like  SPARC)  highlight  a  number  of  challenges  with  implementation  of 
models  for  predicting  properties  over  a  wide  range  of  chemical  structures. 

2.1  Introduction 

This  application  of  green  chemistry  principals  [56]  to  the  design  and  selection  of  new 
energetic  compounds  requires  data  on  their  chemical  fate  properties,  often  in  the  early 
stages  of  evaluation  (including  prior  to  synthesis)  where  experimentation  with  the 
candidate  compounds  is  impractical  [29,  57].  To  overcome  this  obstacle,  robust  models 
are  needed  to  predict  key  environmental  fate  properties,  including  rate  constants  for  nitro 
reduction,  which  is  a  major  detenninant  of  their  fate  [55,  59].  The  reliance  of  these 
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models  on  experimental  input  data  must  be  limited,  due  to  the  challenges  associated  with 
performing  experiments  with  novel  and/or  energetic  compounds.  For  this  reason,  models 
such  as  linear  free-energy  relationships  (LFERs)  and  quantitative  structure-activity 
relationships  (QSARs)  utilizing  calculated  descriptor  variables  (e.g.,  determined  using 
molecular  modeling)  are  desirable. 

2.2  Methods 

Batch  Experiments.  Batch  experiments  were  perfonned  in  60-mL  clear  glass  reaction 
vials  capped  with  butyl  rubber  septa  (Fisher  Scientific).  Each  vial  was  prepared  in  an 
anoxic  glove  box  (5%  H2  in  N2)  and  initially  contained  49  mL  phosphate  buffer,  1  mL 
aqueous  cysteine  solution,  and  an  aliquot  (25-400  pL)  of  aqueous  porphyrin  stock 
solution.  Vials  were  equilibrated  to  25.0  °C  in  a  water  bath.  After  temperature 
equilibration,  the  specified  NAC  was  introduced  to  the  reaction  vial  by  injecting  50  pL  of 
a  0. 1  M  stock  solution.  1  mL  aliquots  were  removed  over  time  and  quenched  by  mixing 
with  an  equal  volume  of  oxygen-containing  (i.e.,  not  deoxygenated)  HPLC  grade 
methanol.  Quenched  aliquots  were  analyzed  by  high  performance  liquid  chromatography 
(HPLC)  using  a  Varian  ProStar  210  solvent  delivery  module,  410  autosampler,  and  330 
photodiode  array  detector,  with  a  Platinum  C18  5p  250-mm  x  4.6-mm  column  (Grace). 
For  most  analyses,  the  mobile  phase  consisted  of  1 : 1  DI  water  and  HPLC-grade 
methanol,  the  flow  rate  was  1  mL  min-1,  and  the  detection  wavelength  was  254  mn. 

Computational  Methods.  Values  of  //'nag  were  determined  from  the  free  energy 
difference,  AGrXn,  for  the  one-electron  half  reaction.  These  energy  differences  were 
calculated  using  electronic  structure  calculations  with  the  COSMO  continuum  solvation 
model  [36]  using  the  NWChem  program  suite  [30].  The  electronic  structure  calculations 
were  perfonned  using  density  functional  theory  (DFT)  [60]  with  the  6-3 1 1++G(2d,2p) 
basis  set  [61,  62]  and  the  LDA  [63],  PBE96  [64],  B3LYP  [32,  33],  PBE0  [65],  and  M06- 
2X  [66]  exchange  correlation  functions. 

Additional  methods  details  are  given  in  Bylaska  et  al.  [67]  and  Salter-Blanc  et  al.  [68]. 

2.3  Results 


2.3.1  Formulation  of  reduction  reaction  pathways/mechanisms  and  mining 
previously-published  data  on  the  reactants  and  products  (Task  2/3.1) 

In  order  to  develop  models  for  predicting  rates  of  NAC  reduction,  it  is  useful  to  consider 
the  mechanism  of  the  reduction  reaction.  In  a  substantial  amount  of  previous  work  (e.g., 
[11,  12,  58,  69]),  LFERs  relating  NAC  reduction  rate  constants  in  the  presence  of  a 
variety  of  reductants  (kred)  to  one-electron  reduction  potentials  (C'nac)  according  to  the 
following  model 


^red 


//I  pi 

=  a  NAC  = 

23RT/F  0.059 


(6) 


were  interpreted  as  evidence  that  the  rate  of  NAC  reduction  is  controlled  by  the  rate  of 
the  initial  electron  transfer  in  an  outer-sphere  reaction  (Equation  7). 
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More  recently,  evidence  from  compound-specific  nitrogen  isotope  analysis  of  NAC 
reduction  reactions  has  suggested  that  the  rate  limiting  step  in  NAC  reduction  is  either  (z) 
the  dehydration  of  y4r-N(OH)2  to  form  Az'-NO,  or  (z'z)  proton  or  electron  transfer  to  the 
nitroaromatic  radical  anion  (yfi -NO2’-),  depending  on  the  reductant  and  solution 
chemistry  [13,  70,  71].  These  conclusions  are  broadly  consistent  with  electrochemical 
studies  of  NAC  reduction,  which  suggest  that  either  dehydration  or  electron  transfer  to 
y4z--NC>2’-  is  rate  limiting,  depending  on  experimental  conditions  [72,  73].  In  order  to 
provide  some  reconciliation  between  the  interpretation  of  LFERs  involving  NAC 
reduction  and  the  more  direct  evidence  from  isotope  fractionation  and  electrochemical 
studies,  and  to  provide  a  more  robust  basis  for  defining  predictive  models,  we  revisited 
the  interpretation  of  the  correlation  between  log(feed)  and  C'nac  as  evidence  that  initial 
outer-sphere  electron  transfer  is  rate  limiting.  This  was  done  using  both  the  model 
described  in  Equation  6  and  a  (non-linear)  free  energy  relationship  (FER)  based  on  the 
Marcus  theory  of  outer-sphere  electron  transfer  [74],  described  for  application  to  organic 
reactions  by  Eberson  [75]. 

2.3.2  Calculation  of  target/response  variable  data  for  training  set  compounds 
with  high-level  theory  (Task  2/3.2) 

Assessment  of  free-energy  relationships  (FER).  In  order  to  reevaluate  and  assess  the 
use  of  FERs  to  describe  NAC  reduction  rates,  and  to  evaluate  the  consistency  of  these 
FERs  with  evidence  regarding  the  rate-limiting  step  of  the  reactions,  we  compiled  the 
major  datasets  for  NAC  reduction  by  homogeneous  reductants  (consisting  of  second 
order  rate  constants,  feed)  and  plotted  the  data  vs.  C'xac  /  0.059  V  according  to  Equation  6 
(in  the  convention  of  previous  studies)(e.g.,  [//,  12,  58,  69])  and  vs.  AG°'  (in  a  Marcus- 
Eberson  plot).  The  results  are  shown  in  Figure  14  A  and  B,  respectively.  The  compiled 
data  include  feed  for  NAC  reduction  by  the  electron  carrier/donor  pairs  juglone/hydrogen 
sulfide  [//,  12],  lawsone/hydrogen  sulfide  [11],  and  Fe(II)  porphyrin/cysteine  [11],  as 
well  as  electrochemically-reduced  9,10-anthraquinone-2,6-disolfonate  (AQDS)  and  Fe(II) 
complexed  by  the  organic  ligands  tiron  [76]  and  desferrioxamine  B  (DFOB)  [77]. 
Throughout  this  work,  we  refer  to  these  redox  systems  simply  as  juglone,  lawsone, 
AQDS,  Fe(II)  porphyrin,  Fe(II)-tiron,  and  Fe(II)-DFOB. 
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Figure  14.  logged)  for  various  NACs  vs.  (A)  £V c/0.059  V  and  (B)  A G°'.  Reductants 
include  the  electron  carrier/donor  pairs  juglone/hydrogen  sulfide  [77,  72], 
lawsone/hydrogen  sulfide  [77],  Fe(ll)  porphyrin/cysteine  [77],  electrochemically  reduced 
AQDS  [70],  and  the  Fe(ll)-ligand  complexes  Fe(ll)-tiron  [76]  and  Fe(ll)-DFOB  [77]. 
Closed  symbols  represent  data  for  which  the  descriptor  variable  was  determined  from 
measured  E\e d.  Open  symbols  represent  data  for  which  the  descriptor  variable  was 
determined  from  E\e a  values  that  were  estimated  from  LFERs  [77,  72],  Only  data  for 
which  E\e d  was  measured  were  included  in  the  fits.  The  correlated  data  were  fit  to  (A) 
Equation  6  or  (B)  the  Marcus  equation.  In  (B),  an  extrapolation  of  the  fit  is  shown  (as  a 
dashed  line)  for  cases  where  the  data  is  well  fit  to  the  Marcus  model. 

The  A1  nac  data  used  in  the  preparation  of  Figure  14  were  those  reported  in  conjunction 
with  the  original  kmi  datasets  [11,  12,  76,  77].  Note  that  in  some  cases,  the  values  of 
Ex nac  were  not  measured,  but  instead  were  estimated  based  on  LFERs  for  juglone  and 
lawsone  [11,  12].  We  included  these  data  in  Figure  14  but  excluded  them  from 
subsequent  fitting.  For  Figure  14B,  A G°'  was  calculated  for  all  NAC/reductant  pairs.  For 
this,  it  was  detennined  that  the  electrostatic  tenn  differentiating  A G°'  from  A G° 
(described  in  detail  elsewhere)  [75,  78,  79]  was  equal  to  zero  in  cases  where  the  reductant 
carries  a  charge  of -1  [75],  and  negligible  in  all  other  cases  [78,  79].  Therefore,  AG°'  was 
taken  to  be  equal  to  AG°.  A G°  was  calculated  from  G'nac  and  the  one-electron  reduction 
potential  for  the  reductant  (E'red). 

2.3.3  Calculate  descriptor  variable  data  for  training  set  compounds  using  lower 
level  theory  (Task  2/3.3) 

Measured  data  for  E]  nac  are  often  lacking  or  difficult  to  obtain,  particularly  for  novel 
energetic  NACs  (e.g.,  due  to  lack  of  synthesized  material  or  experimental  hazards).  To 
overcome  this,  it  is  desirable  to  use  calculated  values,  such  as  those  obtained  using 
molecular  modeling.  Multiple  efforts  have  been  made  to  detennine  A 'nac  using  various 
computational  techniques  [67,  80-82].  We  compiled  and  compared  the  results  of  these 
studies  previously,  in  addition  to  presenting  preliminary  results  from  our  own 
calculations  of  A 'nac  [67].  In  this  work,  we  expanded  on  our  previous  results  to  include 
calculations  using  other  exchange  correlation  potentials  and  applied  these  methods  to  an 
expanded  set  of  compounds. 

A 'nac  were  determined  at  five  levels  of  theory,  using  the  6-3 1 1++G(2d,2p)  basis  set 
with  the  LDA,  PBE,  B3LYP,  PBEO,  and  M06-2X  exchange  correlation  functionals  with 
solvation  detennined  using  the  COSMO  solvation  model.  In  order  to  determine  which 
model  chemistry  gives  results  that  are  most  consistent  with  experimental  data,  a  subset  of 
the  calculated  A 'nac  values  were  compared  to  a  previously  compiled  set  of  measured 
A1  nac  [55].  The  measured  validation  dataset  is  given  in  Supporting  Information  of  Salter- 
Blanc  et  al.  [65],  along  with  the  conesponding  computed  values  detennined  at  the  levels 
of  theory  used  in  this  study.  Deviations  between  the  calculated  and  measured  datasets 
were  characterized  by  calculating  the  mean  absolute  deviation  (MAD),  root  mean  square 
deviation  (RMSD),  and  the  largest  positive  and  negative  deviation  between  the  datasets. 
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Table  4.  Calculated  one-electron  reduction  potentials  (E\  Ac)  and  comparison  to  a  set  of 
measured  values.  Calculated  values  were  determined  with  density  functional  theory 
(DFT)  using  the  the  6-31 1++G(2d,2p)  basis  set  with  the  specified  exchange  correlation 
functionals.  Solvation  was  accounted  for  using  the  COSMO  model. 


Cmpd 

Measured 

Calculated  E1 

NAC  (V) 

Abbr a 

Suatoni b 

LDA 

PBE 

B3LYP 

PBE0 

M06-2X 

B3LYP* 

M06-2X* 

NB 

-0.486 

-0.512 

-0.640 

-0.479 

-0.578 

-0.476 

-0.435 

-0.423 

2-CH3-NB 

-0.590 

-0.659 

-0.771 

-0.616 

-0.695 

-0.630 

-0.506 

-0.528 

3-CH3-NB 

-0.475 

-0.561 

-0.690 

-0.522 

-0.602 

-0.490 

-0.457 

-0.433 

4-CH3-NB 

-0.500 

-0.594 

-0.726 

-0.569 

-0.637 

-0.552 

-0.482 

-0.475 

3-CI-NB 

-0.405 

-0.272 

-0.545 

-0.400 

-0.487 

-0.571 

-0.394 

-0.488 

4-CI-NB 

-0.450 

-0.475 

-0.599 

-0.444 

-0.539 

-0.440 

-0.417 

-0.399 

4-NH2-NB 

-0.568 

-0.964 

-1.043 

-0.753 

-0.863 

-0.737 

-0.576 

-0.600 

3-COCH3-NB 

-0.437 

-0.439 

-0.586 

-0.469 

-0.538 

-0.460 

-0.430 

-0.412 

4-COCH3-NB 

-0.356 

-0.174 

-0.354 

-0.304 

-0.375 

-0.352 

-0.345 

-0.339 

1,2-DNB 

-0.287 

0.033 

-0.169 

-0.121 

-0.230 

-0.299 

-0.251 

-0.303 

1,3-DNB 

-0.345 

-0.339 

-0.581 

-0.390 

-0.505 

-0.384 

-0.389 

-0.361 

1,4-DNB 

-0.257 

0.197 

0.011 

0.059 

-0.075 

-0.117 

-0.158 

-0.179 

2,4-DNT 

-0.397 

-0.542 

-0.638 

-0.486 

-0.618 

-0.404 

-0.439 

-0.374 

2,6-DNT 

-0.402 

-0.318 

-0.718 

-0.569 

-0.682 

-0.525 

-0.482 

-0.456 

TNT 

-0.253 

-0.232 

-0.384 

-0.223 

-0.378 

-0.252 

-0.303 

-0.271 

2-CHO 

-0.355 

-0.151 

-0.345 

-0.271 

-0.361 

-0.341 

-0.328 

-0.332 

4-CHO 

-0.322 

-0.065 

-0.255 

-0.216 

-0.273 

-0.364 

-0.300 

-0.347 

4-CH2OH 

-0.478 

-0.554 

-0.642 

-0.501 

-0.624 

-0.561 

-0.447 

-0.481 

2-ADNT 

-0.417 

-0.586 

-0.729 

-0.554 

-0.642 

-0.466 

-0.474 

-0.416 

4-ADNT 

-0.449 

-0.707 

-0.753 

-0.629 

-0.700 

-0.571 

-0.512 

-0.487 

2,4-DANT 

-0.502 

-0.870 

-0.995 

-0.810 

-0.876 

-0.775 

-0.606 

-0.626 

2-CI-NB 

— 

-0.458 

-0.599 

-0.462 

-0.533 

-0.473 

-0.426 

-0.421 

2-COCH3-NB 

— 

-0.424 

-0.588 

-0.425 

-0.600 

-0.425 

-0.407 

-0.389 

DNAN 

___ 

-0.517 

-0.827 

-0.463 

-0.652 

-0.432 

-0.427 

-0.393 

MAD  c 

0.161 

0.207 

0.099 

0.149 

0.066 

0.043 

0.039 

RMSD  d 

0.043 

0.058 

0.018 

0.031 

0.009 

0.003 

0.002 

Largest  Positive  Deviation 

0.454  e 

0.268  e 

0.316  e 

0.1 82  e 

0.1 40  e 

0.099  e 

0.078  e 

Largest  Negative  Deviation 

-0.396f 

-0.4939 

-0.3089 

-0.3749 

-0.2739 

-0.104  9 

-0.1129 

a)  Abbreviations: 

nitrobenzene 

(NB),  2-methylnitrobenzene 

(2-CH3-NB), 

3-methylnitrobenzene 

(3-CH3-NB), 

methylnitrobenzene  (4-CH3-NB),  3-chloronitrobenzene  (3-CI-NB),  4-chloronitrobenzene  (4-CI-NB),  4-nitroanline  (4- 
NH2-NB),  3-acetylnitrobenzene  (3-COCH3-NB),  4-acetylnitrobenzene  (4-COCH3-NB),  1,2-dinitrobenzene  (1,2-DNB),  1,3- 
dinitrobenzene  (1,3-DNB),  1,4-dinitrobenzene  (1,4-DNB),  2,4-dinitrotoluene  (2,4-DNT),  2,6-dinitrotoluene  (2,6-DNT), 
2,4,6-trinitrotoluene  (TNT),  2-nitrobenzaldehyde  (2-CHO),  4-nitrobenzaldehyde  (4-CHO),  4-nitrobenzyl  alcohol  (4- 
CH2OH),  2-amino-4,6-dinitrotoluene  (2-ADNT),  4-amino-2,5-dinitrotoluene  (4-ADNT),  2,4-diamino-6-nitrotoluene  (2,4- 
DANT).  b)  Dataset  previously  compiled  by  Phillips  et  al.  [83]  c)  Mean  absolute  deviation,  d)  Root  mean  square 
deviation,  e)  Deviant  value  for  1 ,4-DNB.  f)  Deviant  value  for  4-NH2-NB.  g)  Deviant  value  for  2,4-DANT. 
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2.3.4  Task  2/3.4.  Correlation  analysis  and  fitting  QSARs  (Task  2/3.4) 

In  order  to  test  the  use  of  E'nac  values  calculated  using  the  B3LYP*  and  M06-2X* 
methods  as  descriptor  variables  for  QSARs,  previously  reported  values  of  &Fe(ii)P  (also 
shown  in  Figure  14A  and  B)  were  correlated  against  the  newly  calculated  E'nac  values 
(E'nac,  b3lyp*  and  A1  nac,  mo6-2x*).  The  resulting  QSARs  are  shown  in  Figure  15A  and  B 
along  with  linear  fits.  The  strongest  correlation  is  seen  when  using  the  E'nac,  b3lyp* 
values  (r2  =  0.74  vs.  0.55  for  E1  nac,  M06-2X*).  This  is  not  surprising,  since  E'nac,  b3lyp* 
most  closely  approximates  the  measured  values  of  E1 nac  for  the  subset  of  compounds 
used  in  QSAR  calibration  (MAD  =  0.037  V  and  RMSD  =  0.002  V  for  A 'nac,  b3lyp*  vs. 
MAD  =  0.055  V  and  RMSD  =  0.004  V  for  E'nac,  M06-2X*).  The  data  show  unifonn  and 
modest  residuals,  suggesting  that  the  QSAR  calibrated  with  the  A 'nac,  b3lyp*  values  is  an 
adequate  representation  of  the  data.  The  QSAR  resulting  from  the  linear  fit  of  the  data  is 
given  in  Equation  7  (with  standard  deviations  of  the  fitting  coefficients  in  parentheses). 

log  (kmn)p)  =  1 2.4(±2.6)  •  ElNACJ!3L7pt  +  5.8(±1.1) 


Figure  15.  QSARs  for  NAC  reduction  by  Fe(ll)P  calibrated  using  rate  constants 
(log(/cFe(ii)p))  obtained  from  the  literature  [11]  and  E'nac  calculated  at  (A)  the  B3LYP/6- 
311++G(2d,2p)  level  or  (B)  the  M06-2X/6-31 1++G(2d,2p)  level,  both  with  an  applied 
linear  transformation  based  on  correlations  to  a  measured  dataset  (denoted  as  E'nac, 
b3lyp*  and  E'nac,  M06-2X*,  respectively). 

2.3.5  Verification  of  QSARs  with  experimental  values  (Task  2/3.5) 

In  order  to  validate  the  QSAR  given  in  Equation  7  for  predicting  rate  constants  for 
energetic  compounds,  we  measured  rate  constants  for  reduction  of  three  energetic  NACs 
(TNT,  2,4-DNT,  and  DNAN)  by  Fe(II)P/cysteine.  Rate  constants  were  also  measured  for 
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three  non-energetic  NACs  (NB,  4-C1-NB,  and  1,3-DNB)  in  order  to  check  the  agreement 
between  our  results  and  the  previously  published  foeanp  data  we  used  to  calibrate  the  FER 
given  in  Equation  7. 

In  the  absence  of  a  comprehensive  kinetic  model  for  the  cases  that  show  deviation 
from  first-order  behavior,  the  initial  portion  of  the  data  alone  were  fit  to  first-order 
kinetics  (similar  to  the  method  of  initial  rates  [54]).  This  treatment  resulted  in  fitting  one 
to  three  half-lives  of  NAC  disappearance  data.  It  is  likely  that  the  kinetics  at  the 
beginning  of  the  experiment  were  dominated  by  simple  reduction  of  the  NACs,  making 
this  simplification  a  satisfactory  approximation  of  the  rate  constant  for  that  process. 
Details  regarding  data  fitting  for  all  the  NACs  evaluated  are  given  in  the  Supporting 
Information  of  Salter-Blanc  et  al.  [dS].  The  resulting  log(kFe(ii)p)  values  are  reported  in 
Table  5. 
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Table  5.  QSAR  verification  and  validation  data.  /(Fe(ii)p  were  determined  from 
experimental  data. 


NAC 

AfFedDP  (L  M-1  s-1) 

log(#CFep) 

NB 

1.06  ±  0.16 

0.026 

4-CI-NB 

1.946  ±0.083 

0.289 

1,3-DNB 

16.6  ±  0.77 

1.22 

2,4-DNT 

14.34  ±0.76 

1.156 

TNT 

112.1  ±  3.1 

2.050 

DNAN 

17.50  ±0.55 

1.243 

In  order  to  verify  that  our  methods  of  measuring  Areaop  produced  values  similar  to  the 
calibration  dataset,  log(&Fe(inp)  values  for  a  few  non-energetic  compounds  (NB,  4-C1-NB, 
and  1,3-DNB)  were  appended  to  a  plot  of  the  calibration  data  from  Figure  15  (shown  in 
Figure  16)  by  plotting  the  data  vs.  the  associated  E1sac.  b3lyp*  value.  The  data  agree  well 
with  the  QSAR,  thereby  supporting  the  general  agreement  of  our  experimental  data  with 
the  calibration  data.  In  addition,  the  data  for  compounds  overlapping  with  the  calibration 
dataset  (NB  and  4-C1-NB)  show  very  good  agreement  with  the  corresponding  individual 
data  points.  This  suggests  that  the  measured  values  of  /cfC(H)p  are  generally  quite  accurate 
and  that  scatter  observed  in  Figure  16  comes  mostly  from  Zs'nac,  b3lyp*  or  limitations  of 
the  model. 
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Figure  16.  Comparison  of  newly  collected  log(/(Fe(ii)p)  data  to  the  calibration  dataset  and 
QSAR.  Data  for  non-energetic  NACs  were  used  to  verify  agreement  of  the  new  data  with 
the  calibration  data  (method  verification  data).  Data  for  energetic  NACs  were  used  to 
test  the  validity  of  the  QSAR  to  energetic  compounds  (validation  data). 

Following  method  verification,  log(A[  e( u u>)  data  for  the  energetic  compounds  tested  (2,4- 
DNT,  TNT,  and  DNAN)  were  appended  to  Figure  16  to  test  the  validity  of  the  QSAR  to 
energetic  compounds.  The  data  point  for  TNT  plots  in  line  with  the  QSAR  and  the  data 
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points  for  2,4-DNT  and  DNAN  fall  slightly  outside  the  range  of  scatter,  but  is  still  within 
one  log  value.  The  good  agreement  of  the  data  with  the  QSAR  suggests  that  it  is  an 
accurate  model  for  predicting  reduction  rates  of  energetic  NACs  by  Fe(II)P. 

2.4  Conclusions 

In  contrast  to  the  hydrolysis  pathway — where  complexity  in  the  mechanism  and  pathway 
prevented  us  from  obtaining  any  broadly  applicable  QSARs  for  predicting  degradation 
rates — the  effect  of  complexity  in  the  nitro  reduction  pathway  was  effectively  captured 
by  contrasting  two  approaches:  theory-based  correlations  to  El  and  empirical  correlations 
to  fiLUMO.  The  former  provides  a  variety  of  insights  into  the  fundamental  processes 
controlling  nitro  reduction  and  the  latter  provides  simple  QSARs  that  should  be  useful  for 
estimating  relative  reduction  rates  of  nitro  aromatic  compounds.  The  models  were 
calibrated  and/or  tested  with  a  variety  of  model  and  actual  energetic  compounds, 
including  TNT  and  DNAN. 


3.  Objective  4 — Amine  Oxidation 

The  major  results  of  Objective  4  were  published  in: 

Salter-Blanc,  A.  J.,  E.  J.  Bylaska,  M.  A.  Lyon,  S.  Ness,  and  P.  G.  Tratnyek  2016. 

Structure-activity  relationships  for  rates  of  aromatic  amine  oxidation  by  manganese 
dioxide.  Environ.  Sci.  Technol.  [10.1021/acs.est.6b00924] 

Abstract:  New  energetic  compounds  are  designed  to  minimize  their  potential 
environmental  impacts,  which  includes  their  transformation  and  the  fate/effects  of  their 
transformation  products.  The  nitro  groups  of  energetic  compounds  are  readily  reduced  to 
amines  and  the  resulting  aromatic  amines  are  subject  to  oxidation  and  coupling  reactions. 
Manganese  dioxide  (MnCE)  is  a  common  environmental  oxidant  and  model  system  for 
kinetic  studies  of  aromatic  amine  oxidation.  In  this  study,  a  training  set  of  new  and 
previously  reported  kinetic  data  for  oxidation  of  model  and  energetic-derived  aromatic 
amines  was  assembled  and  subjected  to  correlation  analysis  against  descriptor  variables 
that  ranged  from  general  purpose  [Hammett  sigma  constants  (a-),  pKa’s  of  the  amines, 
and  energies  of  the  highest  occupied  molecular  orbital  (Ehomo)]  to  specific  for  the  likely 
rate  limiting  step  [one-electron  oxidation  potentials  (Sox)] .  The  selection  of  calculated 
descriptors  (pKa,  Ehomo,  and  Eox)  was  based  on  validation  with  experimental  data.  Ah  of 
the  correlations  gave  satisfactory  quantitative  structure- activity  relationships  (QSARs), 
but  they  improved  with  the  specificity  of  the  descriptor.  The  scope  of  correlation  analysis 
was  extended  beyond  MnCh  to  include  literature  data  on  aromatic  amine  oxidation  by 
other  environmentally  relevant  oxidants  (ozone,  chlorine  dioxide,  phosphate  and 
carbonate  radicals)  by  correlating  relative  rate  constants  (normalized  to  4-chloroaniline) 
to  Ehomo  (calculated  with  a  modest  level  of  theory). 

3.1  Introduction 

Many  explosives,  pesticides,  and  phannaceuticals  are  nitro  aromatic  compounds,  and 
these  compounds  are  readily  reduced  to  aromatic  amines  under  environmental  conditions 
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[10,  55-55].  Any  source  of  aromatic  amino  compounds  in  the  environment  is  a  significant 
concern  because  this  class  of  compounds  is  characterized  by  high  carcinogenicity  [59, 

90],  Therefore,  the  green  design  of  new  materials  involving  nitro  aromatic  compounds 
must  take  into  account  the  potential  environmental  effects  of  any  aromatic  amines  that 
might  arise  from  their  nitro  reduction.  This  is  challenging  for  the  development  of  new 
energetic  materials,  which  should  be  designed  for  minimal  environmental  impact,  but 
usually  are  not  amenable  to  experimental  measurement  of  environmental  fate  detennining 
properties  (because  of  restrictions  related  to  security  and  safety)  [91-94], 

The  environmental  fate  of  simple  aromatic  amines  has  been  studied  extensively  from 
multiple  points  of  view:  biodegradation,  sorption  to  soil,  wastewater  treatment,  etc.  [95- 
99],  In  most  of  this  work,  the  major  transfonnation  process  is  oxidation,  some  of  which 
leads  to  polymerization  and  oxidative  coupling  to  soil/sediment  organic  matter  [5,  6\. 
Detailed  studies  of  the  oxidation  of  aromatic  amines  are  available  for  manganese  dioxide 
(MnCk)  [100-103],  chlorine  dioxide  (CIO2)  [104],  ozone  (O3)  [105],  alkaline-activated 
persulphate  [106],  iron-activated  hydrogen  peroxide  (Fenton’s  reagent)  [107],  hydrogen 
peroxide  catalyzed  by  polyoxometalate  [108],  triplet  excited-state  organic  sensitizers 
[109,  110],  and  electrochemical  oxidation  [111].  In  this  study,  we  used  MnCF  as  the 
primary  oxidant  system  because  (/)  protocols  are  well  established  for  its  use  in  studies  of 
aromatic  amine  oxidation  [100,  101]  and  (if)  MnCb  is  relevant  to  the  “natural  attenuation” 
of  aromatic  amines  in  soils  and  groundwater,  as  well  as  in  systems  engineered  for 
remediation  [112-116],  A  comparison  of  results  obtained  with  MnCh  to  other  oxidants  of 
aromatic  amines  (CIO2,  O3,  carbonate  radical,  phosphate  radical,  and  triplet  excited-state 
natural  organic  matter)  is  provided  in  the  last  part  of  this  section  on  “cross-correlation” 
analysis. 

3.2  Methods 

Rate  Constants.  Rate  constants  for  aromatic  amine  oxidation  were  measured  in  batch 
experiments  using  MnCF  that  was  freshly  synthesized  using  an  adaptation  of  the  methods 
described  by  Murray  [117]  and  Villalobos  [775].  The  expected  product  from  this 
procedure  is  mainly  bimessite,  which  was  confirmed  by  X-ray  diffraction.  The  resulting 
precipitate  was  washed  three  times  with  deionized  (DI)  water  and  then  suspended  in  pH 
6,  10  mM  carbonate  buffer,  with  0.1  M  ionic  strength  (NaCl).  For  batch  experiments,  the 
MnCh  solutions  were  used  as  prepared  or  diluted  to  the  desired  concentrations  with 
additional  carbonate  buffer,  placed  in  amber  vials,  and  sealed.  Vials  were  then  injected 
with  50  pL  of  10  mM  aromatic  amine  stock  solutions,  shaken,  and  then  stirred  throughout 
the  experiment.  Periodically,  1  mL  aliquots  were  removed  and  the  reaction  was  quenched 
by  removing  MnC>2,  either  by  filtration  or  reaction  with  ascorbic  acid  in  a  methanolic 
solution.  Analysis  of  the  samples  was  performed  by  HPLC  with  C- 1 8  reverse-phase 
column  and  photodiode  array  detection.  Quantification  was  performed  from  absorbance 
peaks  at  234  mn. 

Descriptor  Data.  Several  types  of  descriptor  data  were  calculated  using  methods  that  can 
be  efficiently  applied  to  most  aromatic  amines.  Values  of  pKa  were  calculated  using  the 
pKa  calculators  available  from  ACD/Labs  (as  part  of  Percepta  Predictors  and  SciFinder), 
ChemAxon  (as  part  of  InstantJChem  and  chemicalize.org)  and  ARChem  (as  part  of 
SPARC).  Values  of  F'homo  for  the  aromatic  amines  were  calculated  using  Hartree-Fock 
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(HF)  or  density  functional  theory  (DFT)  using  Gaussian  09  (G09)  [1 1 9]  or  NWChem 
[120]  software.  Free  energies  of  oxidation  were  calculated  for  one-electron  oxidation 
(Eox)  using  NWChem  [120].  For  both  sets  of  quantum  chemical  calculations,  a 
representative  range  of  functionals  and  solvation  models  were  used  and  the  results  were 
compared. 

Additional  details  regarding  the  experimental  and  computational  methods  used  in  this 
work  are  described  in  Salter-Blanc  et  al.  [121] 

3.3  Results 

3.3.1  Formulation  of  amine  oxidation  reaction  pathways/mechanisms  and 

mining  previously-published  data  on  the  reactants  and  products  (Task  4.1) 

Rate  Constants  for  Oxidation  of  Aromatic  Amines.  To  provide  a  training  set  of  kinetic 
data  for  correlation  analysis  and  QSAR  development,  two  types  of  aromatic  amines  were 
selected:  (/)  model  compounds  (usually  anilines  with  one  non-reactive  substituent)  for 
which  kinetic  data  are  available  from  previous  work  on  oxidation  by  MnCF  [100,  101] 
and  (//)  possible  daughter  products  from  nitro  reduction  of  important  energetic 
compounds  [122,  123].  Compounds  of  the  latter  type  are  possible  products  of  nitro 
reduction  at  individual  nitro  group  substituents  on  the  energetic  nitroaromatic  compounds 
2,4-dinitroanisole  (DNAN),  2,4-dinitrotoluene  (2,4-DNT),  2,4,6-trinitrotoluene  (TNT), 
and  3-nitroaniline.  Where  applicable,  the  products  of  nitro  reduction  at  both  favorable 
and  unfavorable  sites  (based  on  computational  predictions  of  regioselectivity)  [1 22]  have 
been  included.  Identifying  information  on  all  of  the  selected  training  set  compounds  is 
given  in  Table  SI  in  the  Supporting  Infonnation  of  [121]. 

In  the  study  of  model  aromatic  amine  oxidation  by  Klausen  et  al.,  [101]  kinetic  data 
were  reported  as  relative  rate  constants  (tad)  using  4-chloroaniline  as  the  reference 
compound  for  normalization.  This  convention  is  analogous  to  normalizing  nitro  reduction 
rate  constants  to  4-chloronitrobenzene,  which  has  been  done  in  many  studies  by 
Schwarzenbach  et  al.  [12,  55,  69,  124]  The  main  rationale  for  this  type  of  normalization 
is  that  it  should  allow  comparison  of  kinetic  data  obtained  in  different  ways  and/or  under 
different  experimental  conditions  (e.g.,  M11O2  mass  concentration  or  molarity,  M11O2  size 
or  surface  area,  and  pH),  thereby  enabling  correlation  analyses  across  more  diverse  data 
sets.  For  the  training  set  of  rate  constants  used  in  this  study,  we  adopted  the  convention  of 
normalizing  to  the  4-chloro  substituted  congener,  and  defined  tael  according  to  Equation 
8, 


(8) 

K  a 

where  tabs  is  the  observed  rate  constant  for  the  compound  of  interest,  and  taci  is  the 
observed  rate  constant  for  4-chloronitrobenzene.  In  Klausen  et  al.,  [101]  only  tael  was 
reported  and  their  data  were  used  here  without  further  adjustment.  Laha  and  Luthy  [100] 
reported  initial  disappearance  rates,  which  we  divided  by  the  initial  aniline  concentration 
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to  obtain  estimated  values  of  tabs,  so  their  data  could  be  used  in  Equation  8.  The  values  of 
tael  for  both  studies  are  summarized  in  Table  6. 


Table  6.  Normalized  rate  constants  (taei)  for  aromatic  amine  oxidation  by  MnCta. 


No. 

IUPAC  Name 

log  kre\ 

Laha  and  Luthy3 

log  kre\ 

Klausen  et  a\.b 

log  kre\ 

This  Studyc 

1 

aniline 

-0.93 

0.48 

-0.100 

2 

3-chloroaniline 

-0.96 

3 

4-chloroaniline 

0.0 

0.0 

0.0 

4 

2-methylaniline 

0.79 

5 

3-methylaniline 

0.79 

6 

4-methylaniline 

0.34 

1.6 

7 

2-methoxyaniline 

1.6 

8 

3-methoxyaniline 

0.68 

9 

4-methoxyaniline 

2.0 

2.5 

10 

3-nitroaniline 

-1.34 

11 

4-nitroaniline 

-4.0 

-  -4.11  d 

12 

2-methyl-5-nitroaniline 

-1.40 

13 

4-methyl-3-nitroaniline 

-1.20 

14 

2-methoxy-5-nitroaniline 

-0.279 

15 

4-amino  benzoic  acid 

0.45 

a)  Calculated  from  kexp  data  reported  in  Laha  and  Luthy  [100] 

b)  Calculated  from  concentration  vs.  time  data  in  Figure  8  of  Klausen  et  al.  [101] 

c)  Details  in  Supporting  Information  of  Salter-Blanc  et  al.  [121]. 

d)  Approximate  value  because  reaction  was  slow. 


3.3.2  Calculation  of  target/response  variable  data  for  training  set  compounds 
with  high-level  theory  (Task  4.2) 

Only  experiment  data  (from  Task  4. 1)  were  used. 

3.3.3  Calculation  of  descriptor  variable  data  with  high-level  and  lower  level 
theory  (Task  4.3) 

Correlations  Based  on  Thomo.  The  highest  occupied  molecular  orbital  energies  (Ehomo) 
of  compounds  often  work  well  as  descriptors  for  correlation  analysis  of  oxidation  rates 
because  oxidation  general  involves  loss  of  electrons  from  the  HOMO  (to  the 
acceptor/oxidant).  Ehomo  is  a  molecular  descriptor  that  is  usually  calculated  from 
chemical  structure  theory  using  commonly  available  software,  but  £homo  is 
approximately  equal  to  the  negative  of  the  ionization  potential  (by  Koopman’s  Theorem 
[725]),  and  the  ionization  potential  (IP)  can  be  measured  by  photoemission  spectroscopy 
and  related  methods. 

For  this  part  of  the  study,  we  were  able  to  compile  experimental  values  of  the 
adiabatic  or  vertical  IP  (VIP)  for  most  of  the  aromatic  amines  in  the  training  set,  and  the 
values  are  tabulated  in  Table  S6  in  the  Supporting  Infonnation  of  [121].  Values  of  Tshomo 
were  calculated  for  all  of  the  training  set  compounds  using  Gaussian09  and  NWChem 
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software;  HF,  B3LYP  or  M06-2X  functionals;  6-3 1G*  or  6-3 1 1++G(2d,2p)  basis  sets; 
and  COSMO,  CPCM,  PCM,  or  SMD  solvation  models.  The  quality  of  the  Ehomo  data 
were  evaluated  by  (/)  plotting  a  matrix  of  correlations  of  all  pairs  of  E\  iomo  data  sets  (not 
shown  because  all  of  the  correlations  were  nearly  1:1)  and  (//)  plotting  the  difference 
between  calculated  Chomo  and  measured  VIP,  as  in  Figure  S2  in  the  Supporting 
Infonnation  of  [ 121] .  Six  representative  sets  of  Chomo  data  were  chosen  for  inclusion  in 
the  figure,  to  reflect  the  full  range  of  calculation  conditions  and  results.  The  most 
significant  dissimilarity  among  the  results  shown  is  that  the  differences  are  relatively 
small  for  Chomo  calculated  with  HF  (<0.5  eV),  but  large  (1-2.5  eV)  for  Chomo’s 
calculated  with  DFT  (M06-2X  and  B3LYP).  This  pattern  has  been  reported  previously 
[126],  and  is  the  reason  that  we  included  HF  calculations  in  this  study.  Within  the  results 
for  each  DFT  method,  the  variability  among  the  differences  are  small  (based  on  MAD, 
RMSD,  and  median),  but  they  follow  consistent  and  expected  trends:  e.g.,  the  largest 
basis  set  (6-31 1++G(2d,2p))  gives  the  smallest  differences. 

Correlations  Based  on  Eox.  The  putative  rate  limiting  step  for  aromatic  amine 
oxidation  is  loss  of  one-electron  to  give  the  corresponding  aryl  amino  radical  [102,  103, 
127],  so  the  free  energy  or  potential  of  this  half-reaction  (ArNH2  ^  ArNH2*+  +  e“)  is  a 
promising  candidate  descriptor  variable  for  mechanistically-based  QSARs.  Data  for  these 
descriptors  can  be  obtained  experimentally — e.g.,  half-wave  potentials  (Em)  from 
electrochemical  measurements — or  by  theoretical  calculations.  Measured  values  of  Em 
from  Suatoni  et  al.  [1 28]  have  been  used  as  descriptors  in  other  QSAR  studies  on  the 
oxidation  kinetics  of  substituted  anilines  [99-101]  (and  phenols  [129-131]),  but  their  data 
only  cover  a  portion  of  the  compounds  included  in  the  training  set  of  this  work. 

Therefore,  we  used  the  measured  values  of  Em  from  Suatoni  et  al.  to  verify  our 
calculated  values  of  one-electron  oxidation  potentials  (Eox)  and  then  used  Eox  for 
correlation  analysis  and  QSARs.  To  obtain  Eox,  we  first  calculated  the  free  energies  of 
one-electron  oxidation  (AG°0X)  with  NWChem  and  10  combinations  of  functionals,  basis 
sets,  and  solvation  models.  From  the  resulting  values  of  AG°0X,  six  representative  sets 
were  selected  and  these  were  then  used  to  calculate  the  final  values  of  Eox  (in  Volts  vs. 
SHE),  which  are  summarized  in  Table  S8  in  the  Supporting  Infonnation  of  [121]  along 
with  the  half-wave  potentials  (£1/2)  from  Suatoni  et  al.  [128]. 

As  we  did  for  the  predicted/calculated  values  of  p Ka  and  E\  iomo,  the  quality  of  the 
calculated  Eox  data  were  evaluated  by  (/)  plotting  a  matrix  of  conelations  of  all  pairs  of 
Eox  data  sets  (again,  not  shown  because  all  of  the  conelations  were  nearly  1:1)  and  (if) 
plotting  the  deviation  between  Eox  and  the  measure  values  of  Em  (Figure  S3  in  the 
Supporting  Information  of  [121]).  The  main  conclusion  from  the  figure  is  that  the  trends 
among  aromatic  amines  are  consistent  between  the  different  levels  of  theory,  and 
therefore  all  should  give  similar  results  when  used  as  descriptor  variables  in  QSARs.  The 
Eox  reported  here  also  correlate  well  (r2  >  0.9,  slopes  ~  1.0,  intercepts  <  0.5  Volts;  details 
not  shown)  with  the  one-electron  oxidation  potentials  for  aromatic  amines  reported 
previously  by  Winget  et  al.  [132],  Arnold  [110],  and  Erickson  et  al.  [109].  Careful 
inspection  of  the  differences  (individual,  median,  MAD,  and  RMSD)  suggests  that  the 
absolute  accuracy  of  E0x  calculated  with  PBE  and  B3LYP  are  slightly  better  than  for 
M062X,  and  that  the  differences  between  COSMO  and  COSMO-SMD  were  negligible. 
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3.3.4  Correlation  analysis  and  fitting  QSARs  (Task  4.4) 

All  6  sets  of  fiHOMO  correlate  well  with  log  kre i  and  the  resulting  QSARs  are  summarized 
in  Table  S9  in  the  Supporting  Information  of  [121],  Consideration  of  the  fitting  statistics 
and  outlier  distribution  in  all  cases  suggests  that  the  6  sets  of  is  homo’s  are  equally 
satisfactory  as  descriptor  variables.  Therefore,  in  the  discussion  that  follows,  we  chose  to 
feature  a  correlation  to  E\ iomo  obtained  using  methods  (B3LYP/6-3 1  G*/SMD)  that  are 
easily  performed  with  common  hardware  and  software  (e.g.,  a  personal  computer  running 
Gaussian09).  The  result,  shown  in  Figure  17,  is  better  than  any  of  the  correlations  to  a 
constants  or  pKa’s,  and  most  of  the  unexplained  variance  (MAD,  RSMD)  is  due  to  just  2- 
3  aromatic  amines  with  the  greatest  residuals.  One  of  these  outliers  is  4-amino  benzoic 
acid  (#15),  which  is  easily  rationalized  because  the  anionic  carboxylate  group  is  likely  to 
influence  the  affinity  of  this  compound  to  the  MnCL  surface.  Another  outlier  is  4- 
nitroaniline  (#11),  which  has  the  greatest  uncertainty  of  the  values  of  kre i  measured  in  this 
study,  due  to  its  slow  reaction  rate  with  Mn02. 


Figure  17.  Correlation  of  log  kre \  for  oxidation  of  aromatic  amines  to  Ehomo.  Ehomo  was 
calculated  using  Gaussian09  with  B3LYP,  6-31 G*,  and  SMD.  Fitting  coefficients  and 
statistics  for  the  regression  line  given  in  the  Supporting  Information  of  [121],  Labels 
correspond  to  numbers  in  Table  6. 

The  relative  success  of  the  correlation  shown  in  Figure  17  is  unsurprising,  given  the 
abundance  of  QSARs  based  on  F'homo  for  oxidation  of  substituted  phenols  (e.g.,  [105, 
133,  134])  and  the  strong  analogy  between  substituted  phenols  and  substituted  anilines  as 
model  systems  for  structure-activity  relationships.  There  are  relatively  few  reported 
QSARs  for  substituted  anilines  based  on  F'homo,  however,  two  are  environmentally 
relevant  and  are  notable  because  they  bracket  the  result  reported  here  with  respect  to  the 
trade-offs  between  rigor  and  scope:  Lee  et  al.  [105]  obtained  more  precise  correlations 
using  rate  constants  for  a  homogeneous  model  system  (aqueous  O3)  and  optimal  selection 
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of  molecular  orbital  energies,  whereas  Colon  et  al.  [99]  reported  a  worse  correlation 
using  rate  constants  measured  in  aerobic  sediment  suspensions  and  Fhomo’s  calculated 
using  semiempirical  methods.  These  correlations  to  Cuomo  are  further  evaluated  in  the 
final  section  on  meta-analysis. 

All  6  sets  of  Cox  correlate  well  with  log  kre  1  and  the  resulting  QSARs  are  summarized 
in  Table  S9  in  the  Supporting  Information  of  [121].  After  consideration  of  the  fitting 
statistics  and  outlier  distribution  in  all  cases  (not  shown),  we  concluded  that  the  best 
correlation  to  Cox  was  obtained  using  the  descriptor  dataset  calculated  using  M06-2X  and 
6-3 1 1++G(2d,2p).  This  correlation  is  shown  in  Figure  18,  along  with  the  regression 
results  that  define  the  QSAR.  Overall,  this  correlation  is  better  than  those  given  above  for 
a  constants,  pKa’s,  or  Chomo’s.  In  this  case,  most  of  the  variance  (MAD,  RSMD)  is  due 
to  just  2  aromatic  amines  with  the  greatest  residuals,  the  overall  confidence  and 
prediction  intervals  are  narrower,  and  the  only  clear  outlier  is  4-aminobenzoic  acid.  As 
noted  above,  the  anomalous  behavior  of  this  compound  is  unsurprising  given  that  its 
carboxylate  group  is  likely  to  alter  the  interactions  of  this  compound  with  the  MnCfi 
surface. 


Eox  NWC/M06-2X/6-3 1 1  ++G(2d,2p)/C0SM0  (V) 

Figure  18.  log  kre i  for  oxidation  of  aromatic  amines  to  E0x.  Calculated  with  NWChem, 
M06-2X,  6-311++G(2d,2p),  and  COSMO.  Fitting  coefficients  and  statistics  for  the 
regression  line  given  in  the  Supporting  Information  of  [121],  Labels  correspond  to 
numbers  in  Table  6. 

The  correlation  between  kK \  and  Em  in  Figure  18  can  be  compared  to  two  previous 
studies  that  used  E0x  as  the  descriptor  for  aromatic  amine  oxidation.  Arnold  [110] 
obtained  correlations  with  similar  statistical  quality,  but  a  significantly  larger  slope,  using 
rate  constants  for  oxidation  of  a  wider  range  of  resonance  stabilized  N-containing 
organics  by  carbonate  radical.  A  rationale  for  this  difference  in  sensitivity  (slope)  is 
introduced  in  the  next  section.  Colon  et  al.  [99]  obtained  statistically  less  satisfactory 
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correlations  for  aromatic  amine  oxidation  using  kinetic  data  obtained  with  a  complex, 
natural  material  as  the  oxidant  (suspensions  of  aerobic  sediment).  Using  MnC>2  as  the 
oxidant,  correlations  that  are  qualitatively  similar  to  Figure  18  were  reported  in  the  two 
previous  studies  and  E1/2  as  the  descriptor  variable  [100,  101],  and  these  correlations  were 
interpreted  as  evidence  that  one-electron  oxidation  was  the  rate  determining  step  in  this 
reaction.  Further  interpretation  of  the  slopes  of  these  correlations  in  terms  of  electron- 
transfer  theory  would  require  converting  E0x  to  free  energies  of  aromatic  amine  oxidation 
by  MnCh  [68],  but  we  were  not  able  to  take  this  step  because  the  necessary  one-electron 
reduction  potential  for  MnCE  is  not  available. 

3.3.5  Verification  of  QSARs  with  experimental  values  (Task  4.5) 

Cross-Correlation  and  Meta- Analysis.  The  approach  to  correlation  analysis  represented 
by  Figure  17  to  Figure  18  is  conventional  in  that  it  describes  the  reactivity  of  a 
homologous  series  of  reactants  (aromatic  amines  with  various  substitutions)  by  one 
reaction  (oxidation  by  MnCfi)  under  similar  or  comparable  conditions  (neutral  pH,  etc.). 
To  detennine  how  these  results  relate  to  other  reactions  or  reaction  conditions,  it  can  be 
informative  to  perform  “cross-correlation”  analysis  between  the  rate  constants  for  one 
reaction  (or  set  of  reaction  conditions)  and  the  rate  constants  of  another  [69,  134].  This 
approach  is  less  common  than  conventional  correlation  analysis,  but  it  can  be  useful  for 
estimation  of  missing  rate  constants  and  for  diagnosis  of  reaction  pathways  or 
mechanisms.  The  best  and  most  relevant  example  of  this  is  the  cross-correlation  analysis 
of  rate  constants  for  oxidation  of  substituted  phenols  by  environmentally  significant 
oxidants  ('02,  O3,  CIO2,  and  H0C1),  which  gives  strong  linear  relationships  for  most 
pairs  of  oxidants  [134,  135]. 

To  extend  this  study  to  include  cross-  and  meta-correlation  analysis  of  aromatic 
amine  oxidation,  we  found  suitable  published  kinetic  data  for  four  oxidants:  O3,  CIO2, 
carbonate  radical,  and  phosphate  radical.  Literature  data  for  triplet  excited  state 
methylene  blue  [109]  and  natural  organic  matter  (3NOM*)  [110]  were  not  included  due  to 
complicating  factors  or  insufficient  number  of  anilines,  respectively.  The  data  we 
compiled,  which  are  summarized  in  Table  S10  in  the  Supporting  Infonnation  of  [121], 
were  subjected  to  cross-correlation  analysis,  but  the  results  (not  shown)  were  of  limited 
value  because  not  enough  compounds  were  common  among  the  data  sets.  To  overcome 
this,  we  perfonned  a  meta-analysis  of  the  data  sets  by  calculating  kK\  (with  Equation  8, 
which  requires  hci,  but  this  value  was  available  in  all  the  major  datasets),  and  comparing 
their  correlations  to  the  best  descriptor  variable  that  can  readily  be  calculated  for  any 
substituted  aromatic  amine  (Ehomo  obtained  with  the  same  conditions  chosen  for  Figure 
17).  The  results  are  shown  in  Figure  19,  with  the  results  from  Figure  17  included,  for 
comparison. 
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Figure  19.  Comparison  of  correlations  for  normalized  oxidation  rates  constants  (kre i)  for 
aromatic  amines  versus  Mn02  (this  study  and  prior  literature)  and  previously  published 
data  for  other  oxidants.  The  data  and  fit  for  Mn02  (gray)  are  from  Figure  17.  The  data  for 
other  oxidants  are  compiled  in  in  the  Supporting  Information  of  [727],  and  the  fit  of  these 
data  is  described  in  the  text. 

The  meta  analysis  represented  by  Figure  19  suggests  three  types  of  trends  in  the 
reactivity  of  substituted  aromatic  amines.  First,  the  data  for  chlorine  dioxide  [1 36,  137 ] 
and  phosphate  radical  [138,  139]  exhibit  no  correlation  to  Ehomo,  which  is  surprising,  but 
the  number  of  amines  are  few,  especially  for  the  phosphate  radical  data,  and  the  original 
reports  of  these  data  provide  only  limited  structure -reactivity  analysis.  Second,  the  data 
for  ozone  [1 05,  140]  and  carbonate  radical  [ 141-143 ]  show  a  very  strong  correlation  to 
Ehomo  for  most  amines,  with  a  small  number  of  outliers  that  are  similar  to  those  that  tend 
to  deviate  from  other  correlations  (e.g.,  nitro  substituents,  etc.).  Fitting  the  ozone  and 
carbonate  radical  data  (pooled,  without  outliers)  gives  the  new  regression  line  (slope  1 .54 
±  0.1 1,  intercept  =  8.57  ±  0.64,  r2  =  0.85)  and  narrow  confidence  intervals  shown  in 
Figure  19.  For  comparison,  data  for  MnCF  (compiled  from  this  and  previous  studies) 
confonn  satisfactorily  to  a  correlation  based  on  .Ehomo,  with  a  slope  (6.07  ±  0.79)  that 
suggests  a  significantly  (roughly  5 -fold)  greater  sensitivity  to  substituent  effects  than  for 
ozone  and  carbonate  radical.  One  reason  for  this  difference  in  slope  is  likely  to  be  that 
M11O2  is  the  only  heterogeneous  oxidant  included  in  the  comparison,  and  interfacial  steric 
or  electronic  effects  on  surface  complexation  might  compound  the  effect  of  substituents 
on  reactivity.  Another  possibility  is  that  the  difference  in  slope  is  the  effect  of  differences 
in  potential  of  the  oxidants  that  is  predicted  by  electron  transfer  theory,  which  would  be 
analogous  to  the  interpretation  we  gave  to  data  on  nitro  reduction  kinetics  in  the 
predecessor  to  this  study  [68]. 

The  concentration  of  data  near  the  fitted  line  for  ozone  and  carbonate  radical  data  in 
Figure  19  undoubtedly  indicates  these  data  are  of  high  quality,  but  probably  also  reflects 
the  relatively  simple  substitutions  of  the  anilines  included  in  these  data  sets  (mostly 
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single  methyl,  methoxy,  or  halogen  substituents)  and  homogeneity  of  the  oxidant 
solutions.  The  training  set  used  in  our  study  includes  anilines  derived  from  energetic 
compounds — which  often  include  multiple  nitro  groups,  and  therefore  are  more 
challenging  for  correlation  analysis — and  heterogeneous  oxidants.  Heterogeneous 
systems,  in  particular,  tend  to  produce  complex  kinetics,  which  is  true  for  our  results  and 
the  two  previous  studies  using  suspensions  of  MnCF  as  oxidant.  Another  study  done  in 
even  more  heterogeneous  systems  (aerobic  sediment  suspensions)  [99]  also  gave  complex 
kinetics,  and  a  correlation  to  E\  iomo  with  a  relatively  high  degree  of  scatter.  For  these 
reasons,  and  because  there  were  few  substituted  anilines  in  common  with  the  other 
studies  included  in  this  meta-analysis,  we  have  not  included  them  in  Figure  19.  However, 
qualitative  analysis  of  the  data  is  consistent  with  relatively  shallow  slope,  similar  to  the 
correlation  for  ozone  and  carbonate  radical,  and  significantly  less  than  the  slope  of  the 
correlation  for  MnCF. 

3.4  Conclusions 

The  ultimate  goal  of  this  meta-analysis  is  to  determine  if  there  is  sufficient  generality  to 
QSARs  derived  for  individual,  model  oxidants  that  useful  predictions  can  be  made  about 
the  reactivity  of  nitro  reduction  products  of  new  energetic  materials  under  environmental 
conditions.  Clearly,  predicting  absolute  rates  of  oxidation  is  not  within  the  scope  this 
work,  but  relative  rates  (e.g.,  as  represented  by  kK i)  vary  with  aniline  substitution  in  a 
qualitatively  similar  way  for  most  oxidants.  These  trends  predict  that  oxidation  (and 
therefore  oxidative  coupling  and  immobilization)  of  putative  reduction  products  of 
energetic  compounds  will  be  relatively  slow  (one-order  magnitude  or  more)  compared 
with  aniline  and  most  common  model  substituted  anilines. 

CONCLUSIONS  AND  IMPLICATIONS 

An  overarching  objective  of  this  project  was  to  demonstrate  that  much  of  the  QSAR 
calibration  process  could  be  done  with  data  obtained  from  computational  chemistry; 
hence  the  project  title  “fully  in  silico  calibration”.  This  could  be  significant  because  some 
emerging  energetics  (e.g.,  insensitive  munition  compounds  (IMCs))  are  not  readily 
available  for  experimental  determination  of  environmental  fate  determining  properties, 
and  therefore  calibration  data  derived  in  silico  can  be  the  only  option. 

With  respect  to  hydrolysis,  the  first  step  in  the  QSAR  development  process  (defining 
the  reaction  that  determines  the  rate  constant)  turned  out  to  be  a  greater  challenge  than 
expected.  All  prior  experimental  and  theoretical  work  on  degradation  of  nitro  aromatic 
compounds  by  hydrolysis  was  compiled  and  shown  to  be  inconclusive  with  respect  to  the 
reaction  pathways  and  products.  New  experimental  and  theoretical  work  perfonned  as 
part  of  this  project  clarified  the  most  favorable  hydrolysis  pathways,  but  there  still  was 
too  much  complexity  to  this  reaction  for  QSARs  to  reliably  describe  the  reaction  rates 
over  a  significant  variety  of  IMCs. 

The  reduction  of  nitro  groups  is  one  of  the  most  widely  studied  environmental 
contaminant  transformation  pathways  and  the  most  important  to  energetic  compounds, 
including  IMCs.  Kinetic  data  on  nitro  reduction  were  compiled  and  analyzed  using  a 
theoretically  based  QSAR  model,  which  allows  one  model  to  describe  rates  of  nitro 
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reduction  by  a  variety  of  environmental  reductants.  Using  new  experimental  data,  and 
newly  calculated  reduction  potentials  of  IMCs  and  related  nitro  compounds,  a  QSAR  was 
developed  that  can  be  used  to  predict  nitro  reduction  rates  of  other  IMCs. 

Using  the  QSAR  for  nitro  reduction,  relative  rate  constants  were  estimated  for  a  large 
number  of  energetic  nitro  compounds,  and  then  representative  IMCs  were  superimposed, 
for  comparison.  Figure  1  shows  that  many  IMCs  will  have  relatively  slow  rates  of  nitro 
reduction,  in  some  cases  more  than  an  order  of  magnitude  slower  than  TNT.  This 
comparison  illustrates  how  QSAR  models  could  be  used  to  guide  the  design  of  IMCs  for 
minimal  environmental  impact. 
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Figure  20.  Correlation  between  the  rates  of  nitro  reduction  for  various  energetics 
compounds  ( k ),  relative  to  TNT,  versus  the  descriptor  variable  that  provides  the  best 
overall  QSAR  (E1).  The  right  axis  shows  rate  constants  converted  to  approximate  half- 
lives  to  give.  The  blue  symbols  are  experimental  data,  obtained  with  IMC  related  amino 
compounds. 


With  respect  to  oxidation  of  amines,  previously  reported  kinetic  data  were  compiled 
and  combined  with  new  experimentally  measured  rate  constants  using  several  IMC 
related  aromatic  amines.  To  develop  QSARs  with  these  data,  many  possible  descriptors 
were  evaluated,  ranging  from  simple  to  advanced.  The  simple  descriptors  gave  adequate 
correlations,  but  the  advanced  descriptors  gave  the  best  QSARs.  The  best  QSAR  was 
used  in  predictive  mode  to  compare  the  relative  sensitivity  of  different  aromatic  amines 
to  oxidation.  In  this  analysis,  most  of  the  IMCs  relate  amines  were  relatively  more 
susceptible  to  oxidation,  which  might  mean  that  these  reduction  products  will  bind  more 
readily  to  soils,  and  therefore  be  less  mobile. 
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